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Using the "Liouville space" (the space of operators) of the massive Ising model of quantum 
field theory, there is a natural definition for form factors in any mixed state. These generalize 
the usual form factors, and are building blocks for mixed-state correlation functions. We study 
the cases of mixed states that are diagonal in the asymptotic particle basis, and obtain exact 
expressions for all mixed-state form factors of order and disorder fields. We use novel techniques 
based on deriving and solving a system of nonlinear functional differential equations. We then 
write down the full form factor expansion for mixed-state correlation functions of these fields. 
Under weak analytic conditions on the eigenvalues of the density matrix, this is an exact large- 
distance expansion. The form factors agree with the known finite-temperature form factors 
when the mixed state is specialized to a thermal Gibbs ensemble. Our results can be used 
to analyze correlation functions in generalized Gibbs ensembles (which occur after quantum 
quenches). They can also be used for non-equilibrium steady states where a constant energy 
flow exists, using the recently derived density matrix. We observe that non-equilibrium form 
factors have branch cuts in rapidity space, we verify that this is in agreement with a non- 
equilibrium generalization of the KMS relations, and we show that the leading large-distance 
behavior of order and disorder non-equilibrium correlation functions contains oscillations in the 
log of the distance between the fields. 
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1 Introduction 

The exact evaluation of correlation functions of local fields is one of the main challenges of 
quantum field theory (QFT). Correlation functions are of crucial importance: they encode, in 
principle, all physical information, and they are immediately related to experimental results in 
condensed matter systems near criticality. 

A particularly interesting family of models of QFT are 1+1-dimensional massive integrable 
models [H El El 0] (see the book [5]). Their extended symmetries allow for exact evaluations 
of many quantities and give rise to a rich mathematical structure to build upon. They have 
also found a multitude of applications to modern condensed matter physics [6] . Much progress 
has been achieved in the past 40 years in evaluating vacuum correlation functions in integrable 
QFT, in particular via the factorized scattering theory and the associated theory of form factor 
expansions [HO [3] (Kallen-Lehmann expansions). The large-distance (or low-energy) structure 
of two-point functions is well understood via the analytic structure of form factors, and both 
their large-distance and short-distance behaviors are under good control (see for instance [TJ, [8] ) , 
giving agreement with expected general principles of QFT. 

In recent years, due to developments in various areas of theoretical physics and to new exper- 
imental results, interest has grown for an accurate evaluation of correlation functions in general 
mixed states. Besides thermal Gibbs states, these also include generalized Gibbs ensembles be- 
lieved to occur after quantum quenches in integrable models [HI EH] (see the review [TT] and the 
analytical results confirming this in the Ising model |12t [T3] ) , non-equilibrium energy-carrying 
quantum steady states [HI [151 ESI EH ) > an d others. In most of the physical applications of cur- 
rent interest, the density matrix is diagonal in the basis of asymptotic states and multiplicative 
on the particles (we will refer to these as diagonal mixed states). There is still relatively little 
known about the structure of two-point functions in such general situations, and how this struc- 
ture depends on the mixed state. Correlation functions in thermal states (or on the cylinder) 
have been widely studied in general QFT (see for instance the books [18] ) , and with more preci- 
sion in massive integrable QFT 0QS[IB[2ll[22l[Ml[M[I3 

133 EH ESI E21 SO]- The structure is relatively well understood, although still much work needs 
to be done in integrable QFT in order to get as powerful large-distance or large-time expansions 
as for vacuum two-point functions. The study of correlation functions in generalized Gibbs 
ensembles is however much more recent, analytic leading- asymptotic results being available in 
the Ising model j!21 |4"T] ; similarly for correlation functions in non-equilibrium energy-carrying 
steady states, there are few results besides a leading-decay result at large distances in the XY 
spin chain |42j . Hence, exact QFT expansions for general diagonal mixed states would shed 
much light onto the structure of correlation functions. 

In the present paper, we obtain for the first time exact, complete large-distance (space-like) 
expansions for real-time correlation functions of order and disorder fields in the Ising model of 
QFT, in general diagonal mixed states. The Ising model is a paradigmatic model of integrable 
QFT, and provides a tool to benchmark many new exact method or new physical context. 
It describes the scaling limit of the anisotropic XY quantum spin chain (including the Ising 
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quantum chain) near the critical value h c of the external transverse magnetic field h. Order and 
disorder fields in the QFT model describe the scaling limit of the order parameter in the ordered 
(h > h c ) and disordered (h < h c ) regimes, respectively, asymptotically closed to the critical 
value. Although the QFT model has trivial asymptotic particle theory - it is a free massive 
relativistic Majorana fermion theory - the order and disorder fields are non-local with respect 
to the asymptotic particles, whence are truly interacting fields. Hence, their two-point function 
should already contain part of the non-trivial structure of two-point functions in (integrable) 
models of interacting particles. Our complete large-distance expansion can be used to analyze 
correlation functions in generalized Gibbs ensembles, where in particular we find agreement with 
the main features of some results of fTl], and also in non-equilibrium steady states, where in 
particular we find agreement with the results of |42j . 

We use the method of form factor expansion with respect to the mixed-state "vacuum" in 
the Liouville space (essentially the space of operators on the Hilbert space). The Liouville space 
construct is an old idea [131 0H 051 061 03 08] which found applications in thermal and non- 
equilibrium physics, and which is ultimately based on the even older GNS construction |49t loTT] 
of C*-algebras (see for instance the book [51] ). However, to our knowledge, it was for the first 
time applied to massive integrable models of QFT in [28]. In the present paper, generalizing 
the thermal situation of |28t [33] , we define the notion of form factors in general diagonal mixed 
states, we calculate them explicitly in the Ising model, and we explain how they can be used to 
obtain large-distance expansions of two-point functions. 

In order to evaluate the usual form factors in the vacuum, the most commonly used method 
in integrable QFT is that by which one solves a Riemann-Hilbert problem in the (analytically 
continued) rapidities of the asymptotic particles [3]. A similar Riemann-Hilbert problem was 
derived in [281133] for thermal form factors in the Ising model, in particular putting in a Liouville- 
space set-up a derivation for form factors on the circle in [27]. In the present case of general 
mixed states, however, these techniques cannot be applied, because we cannot assume any 
strong analytic properties of the density matrix. Further, as explained in [28] (see also [29] ). 
thermal form factors are "analytic continuation" of vacuum form factors in the quantization 
on the circle, and thermal form factor expansions give the real-time version of form factor 
expansions on the circle (compactified imaginary time). Form factors on the circle in the Ising 
model were first evaluated from lattice form factors in free-fermion models on the cylinder 
[521 [53] [Ml E3 ESI EZl [58], and can be evaluated directly in the Ising field theory [261 [27]. Again 
in the present case of general mixed states, these techniques cannot be applied because there is 
(yet) no clear geometric picture behind general density matrices that would lead to alternative 
quantization schemes (but see Subsection 16.21 for ideas concerning a relation with the analytic 
structure of the eigenvalues of the density matrix). 

Instead, our calculation uses novel methods, to our knowledge not used before in the context 
of integrable QFT: we derive and solve a system of bilinear first-order functional differential 
equations for mixed-state form factors. We observe that our mixed-state form factors have 
in general a reduced analyticity region in rapidity space as compared to vacuum form factors 
of integrable models. In particular, in the case of the non-equilibrium energy-carrying steady 
states, the form factors have branch points at zero rapidity. 

It is worth mentioning that, like in [281 133j . our mixed-state form factor method does not 
require any re-summation of partition-function divergencies, involving or not finite- volume regu- 
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larization. These re-summations are required when using the standard vacuum form factors in or- 
der to evaluate, for instance, thermal averages in integrable QFT by explicitly performing traces 
(there is a quite developed technology on this subject [E2[22l[23l[Ml[23[33^ In 
our framework, re-summations are automatically and unambiguously performed in the definition 
of mixed-state form factors (up to their overall normalization), and in the resulting form- factor 
expansion. It is also worth mentioning that our method is a real-time method, something which 
is necessary since as we mentioned, by opposition to thermal Gibbs states, it is not expected 
that there be any clear geometric principle in imaginary time for general diagonal mixed state. 
Our mixed-state form factor method is, in principle, generalizable to interacting integrable QFT 
(although explicit form- factor calculations need more developments). 

Our main observations concerning two-point functions are as follows. Let m be the mass of 
the Majorana fermion, and e~ w ^ be the eigenvalue of the density matrix on the one-particle 
state with rapidity 9 (the normalization is such that the vacuum has eigenvalue 1). If W{6) 
is analytic on the real line, then the equal-time two-point function at separation x > in the 
disordered regime is proportional, in the one-particle approximation, to 

e~ x£ J ^- v (e)e imxsiahe . (1) 

There is a decaying factor where £ is given explicitly in terms of W(8) (see (|55|) ). which gener- 
alizes the thermal-state situation. Further, the amplitude v{9) is: 

e -X-(0) e X+(-$) /-ooiio de , 1 / W(6')\ 

The rapidity contour in (JTJ) can be shifted depending on the analytic properties of W(0). Then, 
the two-point function decays proportionally, up to algebraic or other non-exponential factors, 
to e - x ( £ + msm \ 6 *\ cosh «) as _^ OC) where a + iO* is essentially the singularity of (l — e ±w ^) 1 
making the exponential least decaying. Also, the presence and type of extra factors is determined 
by the type of singularity (pole, branch point, etc.). Similar statements hold for higher number 
of particles, giving the complete expansion, and for any local field. 

If the density matrix represents a non-equilibrium steady state with energy current due to 
left and right asymptotic reservoirs at temperatures 7] and T r respectively, then W{6) is not 
analytic on the real line [14} [T6l I17j . In this case, our analysis suggests that the two-point 
functions have oscillatory terms in log(mrr) with frequencies that are integer multiples of 

1 ( m m \ 

i/ = -]og(orth — tanh— y (3) 

In particular, in the disordered regime, the equal-time two-point function should be proportional, 
at lowest frequency, to 

e -xf 



— cosiv \og(mx) + B) (4) 
mx 

where (in agreement with [42] ) £ ncss = \ (£t; + &r r ) is the average of the leading thermal decays 
at temperatures 2] and T r , and B is some constant. 
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The paper is organized as follows. In Section 2 we review the Ising model (scaling limit and 
field theory), we introduce basic notions concerning the Liouville space and we define mixed- 
state form factors. In Section 3 we present our main results, including the mixing phenomenon 
for normal-ordered product of fermion fields, the exact mixed-state form factors of order and 
disorder fields, and the form factor expansion for mixed-state correlation functions. In Section 
4 we discuss some possible applications of our results: quantum quenches and non-equilibrium 
thermal-flow steady states. In Section 5 we present calculations showing and explaining the 
main results. Finally we conclude in Section 6. Appendices give technical details on various 
aspects. 



2 Ising model and Liouville space 

2.1 Ising and XY quantum chains in the scaling limit 

Consider the infinite-length quantum chain described by the hamiltonian (the generic XY model 

[SUED]) 



" chain ^ / , 



'n u n+l ' 7) u n u n+ 



«+i + ha z n 



(5) 



where a l n are Pauli matrices at site n, k G [— 1, 1] is the anisotropy parameter (in the following 
we will exclude the value k = 0), h is the (dimensionless) external transverse magnetic field 
and J > is an energy coupling constant. The automorphism (a x ,a y ,a z ) t- > (a y , a x , —a z ) 
corresponds to (k, h) i— > (—K,—h), hence it is sufficient to consider the region h > 0. At the 
special values k = ±1 this is the quantum Ising model in a transverse field. This model is 
in general integrable: it is equivalent to a free-fermion model via the standard Jordan- Wigner 
transformation [59l [60], [6Tj . Its phase diagram and critical behavior have been described exactly 

ISa E31 EU - 

The spectrum is described by fermionic excitation modes of energies 



e q = J\/ (1 — n 2 ) cos 2 q — 2h cos q + k 2 + h 2 

for q £ [— 7T, 7r). We see that the model is gapless at the critical value h = 1. This critical point is 
in the Ising universality class for every k ^ 0, hence it is described by the Ising model of conformal 
field theory with central charge c = 1/2 [65] : the free field theory of massless Majorana fermiond3- 
When h approaches the critical value 1, the correlation length £ = — 1| (asymptotically, in 

numbers of sites) diverges. The low-energy excitations are then slowly varying over the chain and 
are correctly described, in the limit, by arbitrary numbers of modes with a relativistic dispersion 
relation (that is, it is consistent to take an arbitrary number of modes and yet take the \ow-q 
form of e q ). One can then construct smooth relativistic fields from which the low-energy and 
large-distance universal behaviors can be predicted. This gives the massive Majorana relativistic 
field theory (see for instance [66]), with Hamiltonian 



H = dx : 



yp Dp - - - 

— it^(x)d x '4>(x) + — i^(x)d x ifj(x) + m iip(x)^(x) : (6) 



1 Often, one only takes a particular locality sector of this theory, which excludes the Majorana fermions them- 
selves. 
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formed out of mutually anti-commutating fields tp(x) and i^{x) satisfying the anti-commutation 
relations 

{ip(x),ip(x')} = {if)(x),ijj(x')} = 5(x - x). (7) 

As usual, in ^ the normal ordering corresponds to the energy of the vacuum state |vac) being 
set to zero. Here vp = \n\Ja is the velocity of the excitations and m = J\h — l\/v F is their mass; 
the constant a is the lattice spacing, a length scale over which ip(x) is expected to be essentially 
constant. 

This field theory will describe correctly physical quantities at energy scales T < mv F = 
J\h — 1| and dimensionful length scales £ > a£ = a\n\/\h — 1|, in the limit where h — > 1 (for 
instance, the energy scale T may be the temperature, and the length scale £ the distance between 
local fields in a correlation function). The full relation is obtained by identifying appropriate 
fields in the Majorana theory with observables in the quantum chain. This identification depends 
on the regime considered: either h < 1 or h > 1, and either k > or n < 0. For simplicity, in 
the following we will restrict ourselves to the region 

k e (o, l] 

and to the observables and cr*. 

These observables are related to three fields in the Majorana theory: two primary twist fields 
a and fi associated to the %i symmetry 

(V,V0^(-V,-V0, (8) 

and the field 

e = i : ifjip : . (9) 

The twist field a, the "order" field, only generates (couples with) even numbers of particles, 
while the twist field /i, the "disorder" field, only generates odd numbers of particles; the former 
has a nonzero expectation value, contrary to the latter. These are fields representing the end- 
point of cuts through which other fields are affected by the symmetry transformation ([5]). There 
are two directions the cut can go on the chain: towards the right, or towards the left. With the 
information of the direction of the cut, we will denote the associated fields by a± and fi±, as in 
[28]. We will provide more detail below. 

The identification in the case of the Ising model, at k = 1, is as follows: 





h-> 1+ 


h->l~ 


a* — (vac | <7^| vac) 


2e(x) 


-2e{x) 


u n 


H+(x) 


a + (x) 



Here we have chosen a direction for the cut; this choice is arbitrary, but must be kept the same 
for all fields in a correlation function. Note that in the low-h regime h — > 1~, <r x takes a nonzero 
expectation value (there is "order"), while in the large-/i regime h — > 1 + , its expectation value is 
zero. In fact, the relation between the chain observables and the quantum fields do not depend 
on the value of k; hence the above holds as well in the general case k £ (0, 1]. The identification 
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implies that, for instance, if n is a quantum-chain observable identified with the field 0{x) of 
scaling dimension d, then, in a state (say translation invariant) of characteristic energy scale T, 

lim {mi) 2d {O n O nl ) T oc (0(x, 0)O(y, 0)) r (11) 

where the limit is taken with x = an, y = an', n — n' = 5£ and T = jmvp for some fixed 
constants 5, 7. Here the QFT fields are evaluated at equal time t = 0, but a similar relation is 
often expected to hold for time-evolved operators. Besides an overall dimensionful factor m 2d , 
the QFT correlation function only depends on dimensionless combinations that do not involve 
the microscopic energy scale J; these turn out to be finite, (x — y)T/vF = 5j and (x-y)mvF = 5: 

(0(x)0(y)) T = m 2d f( 1 ,5). (12) 

The mass m does not vanish in the limit h — > 1 as it is a dimensionful quantity: the limit simply 
indicates that its scale is infinitely far from that of J. In particular, the actual value of m in the 
QFT is here of no meaning. The proportionality constant in is non-universal and depends 
in general on k. Note that the scaling dimension of both a and u is 1/8, and that of e is 1. 

Remark that the quantum chain observable a* — (vac|<r^|vac) cr iti ca i is identified with the field 
±2iip(x)ip(x), similarly to the identification displayed in the first line of (|10p but without normal 
ordering. Its expectation value is negative and logarithmically divergent (this is the logarithmic 
divergence of the shift of the non-critical cr^ expectation value from its critical value). 

In the following, we will set vp = 1. 



2.2 Form factors and correlation functions 

The evaluation of correlation functions of interacting fields in massive relativistic QFT is a rather 
complicated task. In the context of integrable QFT, one may use factorized scattering theory, 
based on the fact that the scattering matrix factorizes into two-particle scattering matrices 
[DElEllI] (see also the book [5]). Then, a very powerful method for vacuum correlation functions 
is that using the spectral decomposition and the exact evaluation of matrix elements involved 
m E2 E]. For a QFT with a spectrum of K particle types, the set of asymptotic states (say in 
states) is described by giving the number of particles N G N in the state, and the associated 
rapidities 6j € R and particle types ej £ { 1, . . . , K } , 

I #1, • • • 5 0iv)ei,...,ejv, 6*i > • • • > On. 

The spectral decomposition follows from the statement that the basis of asymptotic states can 
be used to resolve the identity. From standard relativistic and translation invariance, correlation 
functions in the vacuum at space-like distances r = \J x 2 — t 2 > are then expressed as infinite 
series: 

00 „ 

<vac|0t( M )o(o,O)|vac>= £ £ / d6 l • • • dB N \(y*c\O\0 x , . . . , 9 N ) ei _ eN \ 2 e~ ^ rn S cosh ^ . 

N=0ei,...,e N J0i>->e N 

(13) 

In (|13p . the term with k = is simply |(vac|C|vac)| 2 , and in (vac|0|#i, . . . , 0N)ei,...,e N the 
local field O is implicitly at the space-time point (0,0). This expansion is believed to form a 
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convergent series. Form factors ... ejv (0i, • • • ,@n) m the context of massive integrable models 
of QFT are analytic extensions of matrix elements of local fields (vac|0|0i, . . . ,0N)ei,...,e N an d 
can be evaluated exactly by solving a system of equations and analytic properties, the so-called 
form- factor equations [H Wl\ [3]. We note that form factors are analytic for rapidities in a 
neighbourhood of the real lines, and that the extension from the region 9\ > . . . > On-, where 
they correspond to in states, to regions with different orderings (analytic exchange), gives rise 
to factors of the two-particle scattering matrix; this is Watson's lemma, part of the form-factor 
equations. 

The spectrum of the Ising model contains only one particle type. The space of in states is 
a Fock space T~L over the canonical anti-commutation relations 

{a(9),a(9')} = 0, {a(0), at (<?')} = W ~ 0') (14) 

with the asymptotic states identified as 

|0i,...,0 JV > = at(0i)---a t (0jO|vac) 

(we will denote |0i, . . . ,0jv) = (^li • • • j &n\), where the canonical algebra is in agreement with 
the analytic exchange property of form factors (that is, the two-particle scattering matrix is 
— 1). Naturally, the free fermion fields of the Ising model tp(x,t) and tfj(x,t) can be expressed 
linearly in terms of the operators a(0) and at(0): 

ip(x, t) = J d9 e 9 ' 2 (a(0) ^-iEet + a t (e) e -i Pg x+iE e t^ 

ij(x,t) = J dOe- 6 ' 2 (a(6)e ipsX - iEet -a l! (9)e- ipsx+iEst ^ (15) 



where 



pg = msmh.8, E$ = mcosh0 (16) 



are the relativistic momentum and energy associated to the rapidity 9 and with mass m. The 
Hamiltonian and momentum operators are 

H = J d6E e a ] {9)a{9), P = I d9 p e a 1 \0)a(0). (17) 

The fields a and u in the Ising model are twist fields associated to the Z2 symmetry (|8|). 
Since a couples to even numbers of particles, it is of bosonic statistics; while u is of fermionic 
statistics. With a cut running towards the right we will denote them by cr+ and while with 
a cut running towards the left, by cr_ and Their twist conditions, including the statistics, 
are a direct consequence of the Jordan-Wigner transformation, and are given by 

, \K /\ / ±ip{x')a±(x) (x>x') . . , / ^ip{x')n±{x) {x > x') 

a±{x) ^ x) = \^{x')a ± {x) (*<*') ' ^ X) ^ x) = \±^x')u ± {x) (x<x>) (18) 

and similar conditions hold with the replacement ip 1— > ty. Along with the condition that o± and 
u± be of the lowest scaling dimension, this uniquely defines them up to normalization. 
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Although the Ising model is a free-fermion model, hence has a trivial scattering matrix, the 
order and disorder fields a± and p± are not "free field" , hence have nontrivial form factors for all 
(even or odd, respectively) particle numbers (contrary to the field e). The nonzero form factors 
of the fields a± , fj,± (S3 El EE] and e are given by 



rH^-M (4=) (*) II tanh(^i) 



l<i<j<N 

n^H) = -^sinh(^) (19) 

where (a) := (vac|<r±|vac) is the vacuum expectation value. Note that all other matrix elements 
can be evaluated by using crossing relations [U [6TJ, EJ [8], and that [28] 

al = a ± , nl = ±n±. (20) 

We note finally that form factors of the twist fields a± may equivalently be evaluated using 
Wick's theorem on the particles, with a contraction given by the two-particle form factor. This 
is an indication that a± are normalized exponentials of bilinear expressions in fermion operators 
(the overall normalization is made finite by normal ordering). The twist conditions (|18|) along 
with the requirement that a± be exponentials of bilinear expressions uniquely define, up to 
normalization, the fields a±. Something similar holds for the disorder fields /i±. 

2.3 Liouville space 

Evaluating correlation functions of Ising twist fields in a mixed state is, however, more involved. 
One usually describes mixed-state averages (• • ■)„ via a density matrix p: 

(...) ._ Tr (/°---) 



Tr(p) 



For instance, the (un-normalized) density matrix representing a Gibbs ensemble at a nonzero 
temperature is p = pp := e~@ H , where /3 is the inverse temperature, and with a chemical 
potential there is the extra factor e ~vP f-oc dd ^ ( d ) a ( e \ The density matrix representing a non- 
equilibrium steady state sustaining a constant energy flow is p3] (a similar form has been shown 
in general massice QFT [16], and see also |17] ) 



-A Jo 00 de E g o t(0)a(0)-A, !° x de e b at(e) a (0) 



(211 



where /3 ; _1 and f3~ l are the left- and right-temperatures of the asymptotic baths driving the 
steady state. Further, it has been argued that after quantum quenches, the density matrix 
becomes the exponential of a linear combination of local conserved charges (generalized Gibbs 
ensemble) [TU] , and this has been shown in the Ising model |12l [T5] . Since it is well known 
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that local conserved charges in the Ising model have the form f d9 e sS a)(9)a(9), also in this case 
p is the exponential of an integral over particle densities a'(9)a(9). 

Note that the main elements necessary to obtain the large-distance form factor expansion 
(|13|) are the completeness of the basis of asymptotic states and the fact that the average is 
evaluated in the vacuum. Thanks to the Gelfand-Naimark-Segal (GNS) construction (see for 
instance the book [51]), associated to a density matrix (or more precisely to a state, seen as a 
linear functional on a C*-algebra) there is a vacuum above which one can construct a Hilbert 
space. The resulting Hilbert space is basically the space of operators (more precisely, a certain 
completion of a certain quotient of the C*-algebra), and averages with respect to the density 
matrix are vacuum expectation values in that new Hilbert space. This vacuum expectation value 
can be expanded using the resolution of the identity, so a "form factor expansion" can naturally 
be obtained for mixed-state two-point functions. We use these basic ideas here: we form the space 
of operators on H, which we will referred to as the Liouville space (sometimes referred to as the 
associated Hilbert space) [33l 2H 05] , an d we construct an inner product such that (-) p becomes 
a vacuum expectation value in the Liouville space. This will allow us to obtain a full form factor 
expansion for mixed-state two-point functions. We will not go into the delicate details of how 
to actually construct a C*-algebra in order to mathematically have the ingredients necessary for 
the application of the GNS construction; we will rather provide physical explanations for some 
of the subtleties involved in the process of obtaining a form factor expansion. These simple 
Liouville-space ideas are at the basis of the theory of thermofield dynamics (see for instance 
[4T1 148j). To our knowledge, the first time form factors were considered in the Liouville space is 
in [28} I33j. for the Ising thermal Gibbs state. 

We consider End("H) as the space of operators with basis formed by a ei {9\) ■ ■ -a ejv (0jv) for 
01 > . . . > 9 N , ej = ± and N £ N. Here a+(9) := a){9) and a~(9) := a(9), and we will include 
in this space local fields like cr(x) and p(x), which can be seen as infinite linear combination (we 
omit questions completion and of convergence) . Our Liouville space C p is the inner-product space 
based on End('H), with inner product specified by the density matrix p. With A, B £ End(T^), 
we denote the corresponding Liouville states by \A) P , \B) P respectively, and we set the inner 
product to be 

'<w = (-) 

We restrict ourselves to density matrices p that are diagonal on the asymptotic state basis: 



p = exp 



d9W{9)a\9)a{9) 



(23) 



The function W{9) should ensure that the result is a well-defined density matrix. We will 
consider two cases for it: the untwisted and the twisted cases. Let V(9) be a function of rapidity 
that is integrable on the real line, and uniformly positive: 

inf V{9) > 0. (24) 



In the untwisted case, we simply choose W(9) = V(9). In the twisted case, we consider in p the 
presence of the unitary operator e ln f d9a H s ) a ( s ) that implements the Z2 symmetry ([8]). That is, 

W(9) = V(9) or W(9) = m + V{9). (25) 
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This choice of density matrix includes the usual Gibbs state at finite temperature or chemical 
potential, as well as the non-equilibrium steady state (|21|) and the generalized Gibbs ensembles. 
It may be necessary for convergence of our form factor expansion to also impose that W(9) 
grows fast enough at large \9\, although we will not go into such details (in all examples that we 
consider, it grows at least linearly). 

For convenience, we choose basis elements formed out of the usual annihilation and creation 
operators, but with a particular normalization: 

|vac>" = 1, \e u . . . , 9 N y et _ eN = Q p ei ,..., eN fa, ...,0 N ) a ei fa) ■ ■ ■ <f»{9 N ), (26) 

where the normalization factors are simply related to the Fermi filling fractions, 

N 



QL-,e N (^ ■ ■ ■ > On) ■= II 1 + e-« w M) . (27) 



i=l 

This choice leads to nice analytic properties as will be clear below. In order not to overcount 
basis elements, we need an ordering, for instance 9\ > . . . > 9^. 

Such a structure is sometimes referred to as a "Liouville-Fock" space |46j , and we will refer to 
a doublet (9, e) as representing a "Liouville particle" of rapidity 9 and type e. There is of course 
a hermitian structure inferred from the inner product (the hermitian conjugation still being 
denoted by '). One may evaluate the inner products of basis states using the cylic property of 
the trace and the canonical algebra ([1] 



JV 

J N fa, . . . , 9 N \9[, . . . , 9'^ ^ = J] [(l + e-^WO) g^sfa - flj)] (28) 



i=i 



where we assume the ordering 9\ > ■ ■ ■ > 9^ and 9[ > ■ ■ ■ > 9' N (note in particular that 
the states are not "canonically" normalized). Let a±(0) and its hermitian conjug ate a^(9) be 
operators on C p (sometimes referred to as "superoperators") satisfying the anti-commutation 
relations 

{a e (0),a e ,(0')} = We (9), 4(0} = 0, {^(9), 4,(9')} = (l + e^ w ^) 5^5(9 - 9'). (29) 

The space C p can be identified with the Fock space over this algebra, 

a e (#)|vac>" = 0, {9,,..., fljv)^...,^ = 4 fa) ■ ■ ■ 4 N (9 N )\ vac)". (30) 

The choice of the ordering 9\ > . . . > 9n for describing the basis is important, for instance, 
in the resolution of the identity in terms of basis states. It will be convenient however to define 
\9i, . . . , 9j<[)e ly ..,e N for every 9j € R by the fact that a sign (—1) is obtained upon exchange of two 
Liouville particles (9i,ei) and (9j,ej). Then the resolution of the identity can be re-expressed 
through integrals over the full line, 

°° f°° d9 d9]\[ 

1 = £.5, L mnu\i+e- e > m W 6 ^ ' ' ' 9s)L ~™ *^ {ei ' '" ,dNl (31) 
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Here it is crucial to note that upon such an exchange of the position of Liouville particles, there 
is, by definition, no delta-function "contact" term; by contrast, such contact terms arise upon 
changing the order of operators a + {6) and a~(6') thanks to the anti-commutation relations (|14j) . 
The Liouville space is formed by a continuous basis without discrete (or delta-function) part 
at colliding rapidities, so that the colliding-rapidity submanifold of M> N has zero measure in a 
decomposition of the identity. There is a caviat to this in the case of the twist fields, to which 
we will come back. 

See Appendix [X] for more on the formal structure of C p . 



2.4 Definition of mixed-state form factors 

To every operator A G End(%) one can associate the Liouville left-action 

A £ G End(£ p ) : A e \B) p = \AB) P . (32) 

In particular, the left-action linear map A \- > A £ is an algebra homomorphism, (AB) = A e B e . 
Evidently, averages in the density matrix p are then vacuum expectation values in £ p : 

{A) p = ^(vac|AVac) p . (33) 

Hence, using (j3T|) . two-point functions should have a spectral decomposition on C p where left- 
action matrix elements are involved. This justifies the following definition of mixed-state form 
factors associated to p (generalizing [28} [33] ) : they are (analytic extensions of) matrix elements 

fe P Z,eJ^ ■ ■ ■ M ■■= P (™C\(D%, e N y MN . (34) 

Again, O is implicitly at the space-time point (0,0). These form factors satisfy the relation 

f!';°...,J"i <>.,■<{<■:■■■■ M = -fg°, eri (o 1 ,...,9 j+1 ,e j ,--- ,e N ) (35) 

and the cyclic property of the trace implies 



ei,...,e N 



... , 9 N \0\vacy = f'-:[ , iWv. . . . , (36) 



Mixed-state form factors are essentially traces with insertions of operators a e (6), up to an 
overall factor and up to the subtraction of contact terms at colliding rapidities. Let us make 
this more explicit. Traces of products of creation and annihilation operators are evaluated by 
using Wick's theorem (this follows from cyclicity of the trace and the canonical algebra (|14p ) 
with contractions given by (here with the normalization factor ([2j 



Q £ue2 (6 1 ,6 2 ){a^(e 1 )a^(e 2 )), 



Traces with a further insertion of a local field, expressed themselves through creation and an- 
nihilation operators, are evaluated similarly. This leads in a standard way to a diagrammatic 
expression, associating a single vertex to the local field. Then, form factors are obtained, by 
definition, by summing over connected diagrams, 



f!';°..,,i f h- ...M= few* fa. ■■■^n){0 a ei (^i) • • • a e »(9 N )) p 



connected 
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For instance, 

/f^ 2 (^i^2) = g ei , e2 (ei^2)(Oa ei (^i)a e2 (^))p-(O) p „ (: p 2 ^ 2 |0i) ei . (37) 
Note that the latter can be re-written as 

/ 2 (e 2 \o%y ei = f:f_ e2 (e 1 ,e 2 ) + (0) p f^e^. (38) 

Similar equations generalize (|36p and (|38p to many particles both on the left and on the right. 

A natural physical interpretation of the Liouville space is as the space of particles and 
holes excitations created from the Liouville vacuum consisting of a finite density of particles 
with statistical distribution determined by p. In this sense, a basis element \6\, . . . , 0^)^,...,^ 
represents the presence of N Liouville particles, i.e. particles (e^ = +) or holes (ej = — ) above 
the finite density. A similar construction to that above can in principle be done for general (at 
least integrable) QFT, where these particles and holes can scatter. We do not know yet of a 
clear scattering theory generalizing the usual one to the Liouville space. However, like for the 
usual pure-state case, we may assume that the form factors (|34p . as analytic functions of the 
rapidities continued from the region > . . . > 9n, will allow us to extract this scattering matrix: 
the analytic continuation to different orderings should give rise to particle-exchange properties. 
We will see below that the choice of a minus sign under permutation of Liouville particles in 
the vectors \9\, . . . , 0jv)ei,...,ejv is in agreement with the analytic structure of mixed state form 
factors for the Ising model: the Ising Liouville scattering theory is, expectedly, still that of free 
fermions (but with twice as many particles). 



3 Main results 
3.1 Mixing 

Consider local fields O that are normal-ordered products of fermion fields ip, $ and of their 
derivatives. It is a simple matter to evaluate the traces defining mixed-state form factors (or 
to evaluate any correlation functions) for such fields O by using Wick's theorem as explained 
above. Although these fields are not of immediate physical interest, we believe it is instructive to 
describe a bit more the result of such calculations (see Subsection 15.11 for proofs of the statements 
below) . 

It turns out that the resulting mixed-state form factors of normal-ordered products of 
fermions and their derivatives are simply related to ordinary matrix elements of similar fields 
on the Hilbert space %. For instance, 

(vac|e|0i,0 2 ) (ei = e 2 = +) 



7 2 ,0i|e|vac) (ei = e 2 = -) 
and the zero-particle form factor of the field e (i.e. its expectation value) is 
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(all other form factors are zero). These formulas mean that the mixed-state form factors of ift 
and of e are equal to the matrix elements on T~L of the fields 



lity) := i/>, iX(e):=e + {e) p l (39) 

respectively, where the matrix element taken depends on the configuration of Liouville particle 
types. This phenomenon is in fact completely general. That is, mixed-state form factors of 
normal-ordered products of fermion fields and their derivatives are equal to matrix elements 
on % of linear combinations of such local normal-ordered products with equal and less fermion 
fields. This is a phenomenon of "mixing" , a local field is transformed into a linear combination 
involving its "ascendants" under the free fermion algebra. More formally, there is a mixing map 
il on the space of local fields such that the following relation holds: 

/£,+,_..._ • • . ,0j,e j+1 , ...,e N ) = (e N , . . .,9^(0)^, . . .,e-) (40) 

where there are j plus signs and N — j minus signs as indices on the left-hand side. 

One consequence of ()40p is that mixed-state form factors of normal-ordered product of 
fermions and their derivatives are entire functions of the rapidities, no matter the analytic 
properties of W{9). This is the reason for our choice of the normalization factor (|27j) . 

The map il simply implements the additional "internal" contractions present in evaluating 
the trace that are not present in evaluating usual matrix elements of normal-ordered products. 
It can be represented conveniently using the Liouville space formalism. Indeed, given a field O, 
we can think of the corresponding Liouville state \0) p . Then the map il can be represented by 
the action of an operator on the Liouvile space: 

\n{0))p = v\oy. 

It turns out that the operator U performing the necessary Wick contractions is simply 

a_(0)a+(0)" 



U = exp 



dO 



1 + e w W 



We may also see these additional contractions as occurring as a consequence of a change of 
normal ordering, from one with respect to the usual vacuum |vac) to one with respect to the 
Liouville vacuum |vac) p . If we denote by ° • ° the normal ordering with respect to |vac) p , then 
we have 

\Q)P = °ii(0)*°|vac>'\ 

In fact, all fields il(O) may also be evaluated simply as follows. One first evaluates them for 
all fermion linears and bilinears O; only fermion bilinears get extra contributions, proportional 
to the identity (one internal contraction). Then, the map il on all normal-ordered fields with 
higher number of fermions can be evaluated recursively by using 

mx),li(0)] =il([V(x),0]), mx),ii(0)] =U{$(x) t O]), (41) 

where [•, •] is either the commutator (O bosonic) or the anti-commutator (O fermionic). 
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We would like to mention that the use of the map il in order to evaluate mixed-state form 
factors is in parallel with the exponential conformal change of coordinates in conformal field the- 
ory, which allows finite-temperature correlation functions to be computed from zero-temperature 
correlation functions. But this is here in a massive model, and with general density matrices of 
the form (|23p . It would be interesting to gain a fuller geometric or algebraic understanding of 
it. 

3.2 Exact mixed-state form factors for twist fields 

The map il in principle allows us to calculate mixed-state form factors from the known matrix 
elements on %. However, for the twist fields a and \x, this requires an infinite re-summation: 
twist fields are infinite linear combinations of normal-ordered products (since they have nonzero 
matrix elements for arbitrary large number of particles), hence there are infinitely many internal 
contractions. This re-summation in principle gives rise to two effects: first, the normalization, 
encoded into the expectation value {cr±) p = p (vac|<7j_|vac) p , is modified from its vacuum value 
(a); second, the dependence on the rapidities 9j of the form factors is affected. Here we will 
only consider the dependence on the rapidities. 

Using techniques that do not actually involve the map il, we will show the following (here an 
below, multiple sign possibilities are correlated in any given equation, unless otherwise stated): 

Proposition 3.1 The mixed-state form factors of order and disorder fields are given by 
/f M± (0) 

fc) 

where 

h ^ {9) = \H ^ 9(6 + d0r ' 9{9) = 6XP 
and 

K{6) = -htM- (44) 
i 

Here > is infinitesimally close to ; and the above defines form factors for real rapidities, in 
general, as distributions obtained from boundary values of analytic functions. In particular, form 
factors of o~± and fj,± are analytic functions, except for the kinematic poles at equal rapidities, in 
the strip determined by Im(9j) 6 (0, 7r) for €j = ± and Im(#j) £ (— ir,0) for ej = =F. Further, for 
9 £ M. the functions hf (9) are ordinary integrable functions obtained by continuous continuation 
from these analyticity regions. 

By standard arguments, form factors with higher numbers of particles can be evaluated by using 
Wick's theorem on the particles. The overall normalization is (o~±) p , the contraction of two par- 
ticles ei) and (92, €2) is given by the normalized two-particle form factor fel^{9\,92)/{o-±) p , 
and the remaining single particle (9,e), if any, gives a factor ft^ (9) / (a±) p ; further, there is 



= —ht(9)(a±) p 

= h% { 0,)ht 2 {92) (tanl / 2 - gl± ; (g2 - ei) ° ) £1£2 (a ±)p (42) 



d9 



1 



2™sinh(0-0O 1Og ( tanh ^ 



(43) 
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a minus sign for every crossing of contractions, as is usual for fermions. One can check using 
complex conjugation and (|2(jp that the above expressions are in agreement with (|36|) . 

The functions h^(9) can be referred to as "leg- factors" borrowing a terminology of |271 . The 
fact that mixed-state form factors have the structure of standard form factors times leg-factors 
associated to the individual particles (9j , ej ) is a consequence of the symmetries of the Ising 
model (although we will not make use of this fact in our proof). This was first observed in 
|27| in the case of Ising form factors on the cylinder, but the symmetry arguments are based on 
properties of local fields and on simple dynamical properties of the states, hence they hold as 
well for form factors on the line in any mixed state considered here. That is, the only nontrivial 
part of our proposition is the expression for the leg factors (|43p and (|44D , 

The analytic function h^(9) may be analytically continued from the strip 



Im(0)e(O,7r) (e = ±) 
lm(0) G (-tt,0) (e = T) 



(45) 



where it is analytic, to an extended region, extending on both sides of the strip a Re(#)-dependent 
amount determined by the analyticity of W{6) around the real line. Let us assume that W{9) is 
analytic on some parts of the real line: around these parts, it has regions of analyticity. If 9 lies 
in such a region, and either 9 or 9 ± em lies in If 1 , then the analytic continuation is obtained 
from 

hf(9)hf(9±em) = -^coth^^. (46) 
Further, leg-factors with different values of e are related to each other: 

fc±(0)fc±(0)=±^coth^ (47) 

this being valid for all 9 £ M and all values of 9 in the analyticity region of W{9). This along 
with (j36]) implies that 

h±{0) = ^ih%{d±m) (48) 
whenever the arguments lie in the analytic region of the leg factors. In turn, this implies 

fe P %t.,e N (0i ± eiiMa, ■ ■ ■ M = ±e 1 if p j: i %_ eN (9 1 ,9 2 , . . .,9 N ). (49) 

This can be seen as a kind of "crossing symmetry" . 

Proposition 13. II was derived in [281 [33] i n the case of a thermal Gibbs state, both for W{9) = 
a cosh 9 (untwisted) and W{9) = m + a cosh 9 (twisted) and for any a = m(3 > (where /3 
is the inverse temperature). There, a system of equations and analytic conditions was derived 
using the underlying free-fermion algebra, and its "minimal" solution, assumed to be the correct 
one for the twist fields a± and n±, was shown to be (|42[) (with these W(9)). General QFT 
arguments also suggest a relation between nonzero-temperature form factors and form factors 
on the circle [28]. The latter were independently calculated before \52\ \54\ [27], and the minimal 
solution (|42p was shown to be in agreement [33J. 

However such techniques (analyticity, relation with form factors on the circle) do not work 
for general W{9), because we do not assume any particular analytic properties for W{9) and 
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because there is no alternative quantization scheme reproducing it. Our proof is based on novel 
arguments: we derive from the definitions a system of functional differential equations for form 
factors of twist fields as functionals of W. Its solution is unique once an initial condition, which 
can be taken to be at Re(W(0)) = oo, has been fixed. The form factors given by Proposition 13. II 
do reproduce the correct usual Hilbert space form factors at Re(W(0)) = oo, and we will show 
that they satisfy the system of functional differential equations; see Subsection 15.21 This, in 
particular, provides an alternative proof of the expression for thermal form factors, which does 
not require "minimality" assumptions (but which requires the knowledge of the usual Hilbert 
space form factors). 

In order to gain some intuition into the function (|43p , we will also provide "thermodynamic" 
arguments in order to explain the main features of the leg factors hf(9), see Subsection 15.31 

Analyticity conditions may be used in certain situations beyond the thermal case. We 
will also show that the analytical properties of form factors from Conjecture 13.11 agree with 
a "nonequilibrium KMS relation" that we will derive in Subsection 14.21 for a non-equilibrium 
steady state describing an energy flow; see Subsection 15.41 

Finally, for completeness, we will derive an exact integral-operator expression for mixed-state 
form factors of general exponential of bilinear expressions in Subsection 15.51 



3.3 Mixed-state correlation functions 
3.3.1 Finite expansions 

The resolution of the identity on the Liouville space, expressed in pip , can be used in order to 
obtain a series expression for two-point functions in terms of form factors (|34]) . leading to an 
expansion similar to (|13p . Formally, then, taking into account the state normalization (|28p . we 
expect to have 

( (M) t ( o,o)>, = £ £ / ■ 

jV=0ei,...,6 JV - / - 00 llj=ll i + e 3 

(50) 

Whenever O is a field whose form factors are zero for large enough numbers of particles, for 
instance the fermion fields tp and ip or the field e, then this is indeed a correct expression. In 
these cases, as we noted, form factors are entire functions of the rapidities. The integral over 
9j is convergent if ej < 0; it is conditionally convergent if ej > 0. Assuming without loss of 
generality that x > 0, in the latter case the integral can be made into a convergent integral, for 
space-like distances (x 2 > t 2 ), if W(0) is analytic on parts of the real line far enough from 9 = 0. 
Indeed, this is done by the shift 6j i— > 6j + irj in the region |Re(#j)| > K, for K > large enough 
and r\ > small enough in such a way that 9j remains in the analyticity region of W(9j). 



3.3.2 Twist fields two-point functions 

For twist fields o± and (i±, the form factor expansion is infinite. It turns out that (|50p is modified 
in various ways due to the presence of the semi-infinite cuts that emanate from their position. 

First, consider the general form of the expansion ^ s p (vac\0(x, t)\s) p p (s\O Jf (0, 0)|vac) p . If 
this is to be a large distance expansion, then we should interpret the intermediate states \s) p p (s\ 
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as lying, or carrying correlations, over the region between the fields, and the vacuum states p (vac| 
and \vac) p as representing what is happening far on the right and left, respectively (recall x > 0). 
Where there is a cut due to a twist field, the density matrix is effectively modified, p i->- pr, with 

W{9) h+ W*{0) := W{9) + iir. (51) 

It was argued in [33], in the case of the finite-temperature form factors and via a comparison 
with form factors on the circle, that the intermediate states must lie in a region that is not 
affected by the cuts. This is in order for the space of excited states constructed, and the form 
factors obtained, to correspond to the sector represented by p, and not the modified one pK 
Hence, again in the space-like region x 2 > t 2 and with x > 0, form factor expansions can be 
obtained for 

(a+(x, t)a- (0, 0)) P) (p+(x, t)p-(0, 0)) p ; 

note that the cuts are emanating away from the region between and x. On the other hand, 
the scaling limit of lattice correlation functions of the cr^ Pauli matrices, in a state pi a t that is 
the lattice correspondent to p, naturally lead (see the table (fTU|) ) to the correlation functions 
{cr+(x, i)<7_|_(0, 0)) p (ordered regime) or (p + (x, t)p + (0, 0)} p (disordered regime): in the Jordan- 
Wigner transformation, the Pauli spin matrices are written in terms of infinite products of 
fermion operators starting at the matrix's site and going in a fixed direction. That is, 

,. , A \2d/ x x \ i / n \ / nw f a (ordered regime) , roN 

\\Y&(miY a (.olal,) , . oc (u)+(x,0)u)+(y,0)) o , u = < ) ,. , , 6 . ' (52) 
h^V J s n n /plat v +v 1 y " p \ p (disordered regime). v ; 

Hence in order to obtain form factor expansions for these, we need to use 

t)uj + (0, 0)) p = (u}+(x, i)w_(0, 0))ptt (oj = cf or u = p). (53) 

Note that these correlation functions should be real. 

Second, the presence of a cut affects the translation covariance property of x-dependent 
matrix elements. This property can be expressed as 

"(vaciM*, t)\e u . . . , e N y ei _ eN = < r " ' i: '-'\f !';:...,, ^- •••,%) (54) 

(where u = a or u = p), with 

— m cosh 8 log I coth I . (55) 

- oo 2tt V 2 / 

The factor e^ x£ can be understood physically by the fact that the trace in the region where a 
cut lies contributes a different free energy to that where no cut lies, so that, after appropriate 
normalization with respect to the point x = 0, a shift of the cut end-point gives rise to the 
exponential of the twist-field free energy deficit. We will provide a thermodynamic argument 
for it in Subsection 15.31 Note that if is associated to W", we have 

£" = -£. (56) 
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The final subtlety is that the form factors of twist fields are not entire functions of the 
rapidities - they are distributions defined as boundary values of analytic functions, and it appears 
to be important, in order to obtain a well-defined form factor expansion, to be able to shift 
contours towards the analyticity regions. Hence, we need to further require that W(9) be 
analytic on a neighborhood of R. This requirement will be carefully lifted in the particular case 
of the non-equilibrium steady state (|2"T|) . 

Putting these subtleties together, we obtain the following. 

Conjecture 3.1 With lu = a or u) = fx, and with W{9) analytic for 9 £ R, we have 
{u> + (x,t)uj + (0,0)} p (57) 

N=0 ei,...,ejv ' llj=l I 1 / 

Recall that is the density matrix corresponding to in 1151}) . 

We used (|53p , then from this the expansion obtained from , then (|54p and (|36p , and finally the 
t( modifications given by ()56p and ()51|) . As before, in order to obtain a large-distance expansion 
from the conditionally convergent integrals, we can shift the 9j contour in (|57p by ieji] for 7] > 
small enough in such a way that #j remains in the analyticity region of W(9j). Note that these 
shifts keep the rapidities in the analyticity region of the form factors involved. This means that 

/ , nx fnn \K -x£, \ / \ f l + 0(e- 2mxsin l e *l cosha ) (u = a) rnQ , 

K(x,0)w + (0,O)) p = e (<T + ) ncss (°-)n CSS ■ I 0(e _^sin r |cos ha) L = /i) ^ 

as x — > oo, where z := q + z^*, with 0* G [— 7r/2, 7t/2] \{0}, is the position on the upper boundary 
of the analyticity region of (l- e -wW)-i , or on the lower boundary of the analyticity region of 
(1 — g^w) - ^ making the exponential least decaying (note that sin|#*|cosha = |Im(sinh z)\). 
Here, 0(e~ kmxsm \ e *\ cosha ) includes possible algebraic or other non-exponential factors in mx. 
These factors are determined by the type of singularity at a + i9*. We provide examples of this 
in our analysis of a particular generalized Gibbs ensemble in Sub-section l4.lt and in our analysis 
of the non-equilibrium steady state in Sub-section 14.21 

We may re- write (|57p explicitly, for any sign of x (in the space- like region) , as 



<<7 + (x,t)<7 + (0,0)} p ± </i + (x,t)//+(0,0)> 




n 



u' =1 



j<k 



where 



tanh J ' o v 3 1 (59) 

^-^(e + ieO)" 26 



K e {9) 



2ir (1 - e-^W) 
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and g(9) is given in (|43p . We provide in Appendix[B]a check to one-particle order that the result 
is a real function, as it should be from general arguments. 

Conjecture (|3.ip was shown to be correct in the thermal Gibbs state in |33j by showing, 
through contour shifts, that it reproduces the correct form factor expansion in the quantization 
on the circle. 

As we mentioned, in the case of the non-equilibrium steady state, because of the non- 
analyticity at 9 = 0, care must be taken in performing the shifts, and this affects the large- 
distance behaviour (see Subsection I4.2|) . In particular, although the form of the leading expo- 
nential decay is still valid, the subleading terms are slower than exponential. Note that the factor 
e~ x£ is in agreement with the result of [42] in the anisotropic XY model out of equilibrium. 

4 Applications 

4.1 Quantum quenches 

Quantum quenches are physical setups where, in the usual convention, one suddenly modifies one 
or more parameters (couplings, etc.) of a system, usually initially in its ground state, and let it 
evolve. Hence in general, one has two hamiltonians H and H' which do not commute with each 
other, and the ground state of H (or some low-energy excitation, or some natural mixed state) 
is unitarily evolved with the evolution operator e~ ltH . The main problem is to determine what 
happens after a large evolution time t; for instance: is there thermalization? A ground breaking 
experiment |68j showed that dimensionality plays a crucial role, and suggested that conservation 
laws may have a strong influence on the result. A large body of work has been done in this 
subject (for a comprehensive summary, see for instance |llj). When measuring averages of local 
operators in infinite volume, one would indeed expect the large-time results to be consistent 
with thermalization, with a temperature that depends on the actual quench performed and 
on the initial state. It has been observed that, starting with ground states, thermalization 
indeed occurs in non-integrable systems, but that the infinite number of conserved quantities in 
integrable systems restrict the dynamics enough so that usual thermalization does not occur. 
In integrable systems, the final state is rather one described by a generalized Gibbs ensemble 
(GGE) PHI)]. A GGE has a density matrix of the form 

/>GGE = e- E -i ftlH " (60) 

where H n are the local conserved quantities, and j3 n are the associated generalized inverse tem- 
peratures. The generalized temperatures are determined by the requirement that the averages 
of the conserved densities in the GGE be equal to those in the initial state. Note that we only 
expect conserved quantities that are bounded from below to appear, and usually there are no 
flow of energy, particle, etc. In these cases, although the state is not strictly at equilibrium (be- 
cause we do not have a standard Gibbs' ensemble), it is in a natural "generalized" equilibrium. 
In the Ising model, a precise study of quenches shows that indeed the GGE occurs \12 \ I13 |. [4TTj . 

Our results can be applied directly to GGEs in the Ising model thanks to the trivial observa- 
tion that in the Ising model, conserved quantities are linear combinations (integrals) of a^(9)a(9) 
(a similar statement holds for general integrable QFT). One usually may construct generalized 
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hamiltonians using fixed-spin charges Q n via H n = (Q n + Q_ n )/2, with 

Q n = J dd e nd a\6)a(6) (61) 
(and in particular H\ is the usual hamiltonian) , and one has, according to (|23p . 

oo 

W(0) = ^Pn cosh(n#). (62) 
n=l 

Naturally, more general linear combinations of the Q n s are possible. Hence our results (|5Up (for 
fermion fields and all their descendants) and (|57p (for order and disorder fields) give, in principle, 
the full large-distance expansion of correlation functions in the Ising model in any GGE. 

A particular quench corresponding to an abrupt change of the external magnetic field from 
ho to h was considered in [12[ [T3l I41j . In the scaling limit, with both ho and h near to 1, this 
corresponds to a change of mass from mo to m. The result of these works can be expressed, in 
the scaling limit, by the following choice of W{6): 

tanh ^ = Sinh2g + ??mo/m =: U (63) 
cosh 0a/ sinh 2 9 + m^/m 2 

where rj = + for a quench from ferromagnetic to ferromagnetic or from anti-ferromagnetic to 
anti-ferromagnetic regimes, and r\ = — for the other cases. We see that for r? = — this is out of 
the context that we considered, since then W(9) < for small enough values of 9. However for 
7] = + this can be treated with the present formalism. In particular, we find that 

Both of these have poles at 



and branch points at 



sinh 9 = ±i sj mo/m 



sinhfl = iimo/m. 



Hence, this implies that our expansion (|57|) (or (|59]> ) provides a full large-distance expansion for 
correlation functions of ordered and disordered fields in the universal stationary regime occurring 
after such magnetic-field quenches, and in particular that the leading next-order decay is (|58p 
with 

sin 1 0*| cosh a = \J mo/m or mo/m. 

We may perform a short analysis in order to show that our results agree with the general 
features of the results of [41] . First, in the ordered (ferromagnetic) regime, we must use the 
order field <r+. According to ([58]) . we then see a pure exponential decay controlled by £ . Indeed 
a pure exponential decay was derived in [H]. In the disordered (anti- ferromagnetic) regime, 
we must use the disorder field //+. We now see in ([55]) an exponential decay, with a possible 
algebraic factor. In order to evaluate the possible algebraic factor, we simply need to deform 
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the 9 contours, and look at the singularities, at points z, making |Im(sinhz)| minimum. If 
mo/m > 1, then ^Jmojm < mo/m, hence the poles are these singularities. Poles provide a 
pure exponential behavior, so we have a pure exponential decay. A pure exponential behavior is 
indeed what was found in [3T] for ho > h in the anti- ferromagnetic regime. If mo/m < 1, then 
y^mo/m > mo/m, and the branch point is the nearest singularity. Hence we have an exponential 
decay, but branch cuts give additional algebraically decaying factors. Indeed, an algebraic decay 
is obtained in [U] in this situation as well. The branch cut is of square-root type, so that, after 
our contour shifts, we have to evaluate integrals of the type 

poo 

e~ mox / dp ^//3e~ Cmx oc (mx)- 3/2 e'™ *. 
Jo 

An algebraic decay with exponent —3/2 was indeed found in |41j . 

We may also describe the correction terms, which were not described in [3T]. We should 
have for instance, in the ordered regime, considering the two-particle contributions of the form 
factor expansion, the main exponentially decaying behavior times 

1 + O ((mx)~ 3/2 e -(V^^+^o)A or dered, m > m 
1 + ((mx) -5 e -m ° x ) ordered, mo < m. 

A similar, but simpler, analysis may be done for the field e, using the finite-sum exact 
expression ([50]) . 

We note that our context may in principle be generalized to cases with W{9) < 0; the 
difficulty to overcome will be to deal with singularities in the denominator of (1571) . 

Hence, a simple form factor analysis provides the main features of some results of [5T] as well 
as some simple corrections. A detailed analysis of the large-distance expansion from our exact 
result (|57p . further confirming and generalizing |41j . would be very interesting, but is beyond 
the scope of this paper. 

4.2 Non-equilibrium steady states 

A non-equilibrium steady state (NESS) can be seen as a state which is steady with respect 
to the dynamics of the model, but where there are flows of energy, particles, charge, etc with 
constant rate. Non-equilibrium steady states can be "created" by evolving unitarily, for an 
infinite time, two semi-infinite halves of a system initially separately thermalized. If the two 
halves are thermalized at different temperatures Pf 1 and /3" 1 (left and right halves respectively), 
but no bias in any other chemical potential, then the steady state obtained possesses energy flows 
but no other flows. This state was exactly described in critical systems in [15] using conformal 
field theory, near to criticality in |16] using general massive quantum field theory, and in the 
XY quantum chain in [T3] (see also [IZ])- The latter two results can be stated, in our present 
notation, by saying that the function W{9) is specialized to 

W ness (9) ■= frE e e(8) + (3 r E e e(-8) (64) 

where Q(x) is the Heavyside step function (we denote by p ness the associated density matrix). 
This associates an inverse temperature to right-moving particles, and an inverse temperature 
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/3 r to left-moving ones. Physically, this formula is justified by the fact that right-moving particles 
come from the far left, which was thermalized at temperature ft -1 , and vice versa. 

Note that we may also consider the twisted density matrix p| e ss- Although this does not 
have as clear a physical meaning in the context of the Ising spin chain, it is this twisted density 
matrix that is needed in the form factor expansion (|57p of the physical correlation functions. 



4.2.1 Non-equilibrium form factors have branch cuts 

Since the function W ness (6) is analytic on Re(0) < and on Re(#) > but not at 9 = 0, it 
does not satisfy the assumptions that we made for the form factor expansion. This means that 
we can analytically continue form factors fei, c .^,'e^ ■ ■ ■ ,&n), for U) = (J, or uj = a, from their 
analyticity strip 9j G if. (|45p whenever Re(#) < or Re(#) > using ()46p . but not from the 
point 9 = 0. For /?; ^ {3 r , we can see from (|46p that there is a jump, in this analytic continuation, 
in going from positive to negative real part of 9. This means that the functions hf(9) (hence 
also all form factors) have branch points at 9 = and 9 ± ein with branch cuts running away 
from the analyticity region. The jump through these branch cuts can be evaluated using 
and ([64"]) . and is given by 



hf(9-0) 
ht(9 + 0) 



m/3 r cosh 9 mfii cosh 9 
coth tanh (for VK ncss ) 

(65) 

m/3 r cosh 9 mfc cosh 9 » 
tanh coth (tor W-ness) 



for Re(6>) = and 9 ± em G if. In fact, one can evaluate the leading small-0 behavior by 
extracting the pole 1/(9 — 9') from the factor l/sinh(# — 9') in (|43p . thus obtaining 

h ±, a , / ° Tin (forW^) 1. / rnp r mft\ 

^ W K ( ff**r ( for wLs) ' 7 := ^ l0g l C ° th ~2~ tanh — ) ■ (66) 



4.2.2 Expansion for two-point functions: oscillations in log(mx) 

We can attempt to make the expansion (|57|) into a large-distance expansion by shifting contours. 
This is obtained by shifting the ^-contours towards the positive imaginary direction by an 
amount +iir, for every rapidity associated with ej = +. In shifting the contour, poles are picked 
up from the factors 1 — e~ ejWncas ^ 9j ^ in the denominator; and the segment running from to m 
is surrounded in such a way that the contour passes by 0, because the factors 1 — e~ ejW ^ 9j ^ are 
not continuous at 0. Fully shifted contours with e = + cancel unshifted contours with e = — 
thanks to (|49|) . so that one is left with pole contributions and integrals along the segment from 
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to in. The result is (see Appendix [C] for some details) 



(u+(x,t)cj + (0,0)) Pne 



xSnc 



N,M=O ni ,...,n N eZ\{0} 



n d9\ ■ ■ ■ d0M \ TV .M 



N\M\ 



(2^) JV i M x 



N —mx cosh(a(nj))+2irnjt/ /3(nj) 



n 



L mj3(rij) cosh a(n^ 



M 

n 

fc=i 



-mi sin(6fc)— imt cos(0fe) 



1 



1 



x F ( a(m) + —,... ,a(njv) + yj^i> ■ ■ ■ ,Wn 



where 



and 



a{n) = arcsinh 



27rn 



ft 



(n > 0) 
(n < 0) 



(67) 



(68) 



(69) 



s /3(n)m / 

Here again, either uj = a ov oj = jjL. 

We see that the standard Matsubara frequencies of finite-temperature correlation functions 
appear on the second line of the right-hand side in (|67|) ; however here there are two tempera- 
tures, hence two frequencies, depending on the sign of rij. The Matsubara frequencies can be 
interpreted as coming from the quantization of the momentum in a quantum system on the 
circle, as a finite-temperature state can be interpreted as a condition of (quasi-)periodicity in 
imaginary time, with extent given by the inverse temperature. Here, we see that we have the 
remnant of this in our non-equilibrium steady state, but with two different circumferences, for 
right- and left-moving particles. This does not make full sense in the massive theory, as right- 
and left-movers do not separate (fields do not factorize). This is why we have extra integrals 
on finite intervals; these can be interpreted as providing the bridge for right- and left-moving 
modes to jump from the circle of circumference ft to that of circumference ft and vice versa. 

The leading exponential decay e~ x£ncBB is in agreement with the result of [32] obtained in 
the (anisotropic) XY model out of equilibrium, after one performs the scaling limit. The non- 
equilibrium "free energy deficit" £ ncss can be written as the average of the free energy deficits 
£p associated to equilibrium thermal density matrices, 



ftlCSS - ^ + £,9r) > £/8 : = 



00 dO , „, / , m/3cosh(#) 

— m cosh log coth 

.~ 2?r 6 V 2 



(70) 



The form of the leading or subleading terms multiplying this exponential decay can be 
obtained from the one- or two-particle integral and zero-particle sum, M = 1,2 and N = 0; 
these are not exponential, hence we can neglect all exponentially decaying terms coming from 
higher values of N . They are evaluated by approximating the integral over [0,7r] ofe- mxsin ( 6 h{6) 
by the integral over [0, oo) of e~ mxS (c(9) + c(ir — 9)). Let us take t = for simplicity. Using (|66p 
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and (|46p . and the fact that the two-point function is real, this gives in the case of the disorder 
field (M = 1): 

</i + (x,0)/i + (0,0)> PncBB =e-* e °~ (<J + ) PncBB (*-) PncBB (—cos(2 1 log(mx) + B) + ..) (71) 

where 7 is defined in (|66p . The constants A and B can also be evaluated from the M = 1 
integral after a somewhat tedious analysis (see Appendix O for the answer). However, we do not 
expect this M = 1 evaluation of the constants A and -B to be meaningful. Indeed, the omitted 
part in (|7ip . corresponding to higher values of N and M, contain terms of exactly the same 
form as the first subleading term written in (|71|) . but with different constants A and B, with 
A possibly logarithmically divergent. These come from higher values of M still with N = 0, 
and arise due to the "kinematic poles": the singularities of the coth((6q — 62)/ 2 factors - see 
Appendix [C] for some details. Yet we expect that the oscillating form obtained from the M = 1 
analysis be correct, and the possible divergences to be re-absorbable into the normalization 
of the field. Note that other subleading terms are terms with higher oscillating frequency in 
log(mx) and with higher power of l/(mx) (coming from higher values of M with iV = 0), as 
well as exponentially decaying terms. 

A similar analysis for the order field is obtained from M = 2. In this case, the first subleading 
term is 0(1); this is for the same reason as that for which the next subleading term in the 
disorder-field case are of the same form as the first subleading term. Hence, also in this case, the 
vacuum expectation value (a±) Pncaa does not reproduce the correct normalization for the leading 
exponential decay of the two-point function, (a+(x, 0)cr + (0, 0)) Pncss ~ C e~ x£ncBB with in general 
C 7^ {&+) Pncss (^Opnoss- Here the full physical meaning, beyond its involvement in the form factor 
expansion ([67]) . of the "expectation value" {<J±) PncBB is not entirely clear - such expectation values 
are usually determined from the large-distance asymptotic under the condition of conformal 
normalization at small distances. 

An important conclusion, however, is that in general, both for the order and disordered 
regimes, we expect to have oscillatory terms in log(mx) with frequencies that are multiple of 
27. Yet it would be interesting to understand more fully the large-distance expansion. 

4.2.3 "Non-equilibrium KMS relation" and analytic properties 

In an equilibrium Gibbs state at temperature the two-point function of fermion fields 

Gp(x,r) = (tp(x, r)V>(0, 0)) P/3 at imaginary time t = —it, t G R satisfies the quasi-periodicity 
relation (KMS relation, see for instance [18j ) 

Gp(x,T) = -Gf}(x,T + P). 

This relation is obtained by using the cyclic property of the trace and the relation e^ H ip(x, r)e~^ H = 
tp(x,T + 0). It is at the basis of the description of the Gibbs state in terms of Euclidean 
(imaginary-time) fields on a cylinder of circumference j3 (for fermions: with anti-periodic con- 
ditions). This relation is however broken out of equilibrium. In the case of the present free 
fermion theory there is a straightforward "nonequilibrium KMS relation". Recall the mode ex- 
pansion (|15|) for the fermion fields. We may define non-local fermion operators by integration 
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over rapidity half- lines, 

i/; l ' r (x, t) = - x f^ f d9 e e/2 (a{9) e ipsX - iEet + aUO) e - ipeX+iE <>] 
2 V Je^o v 1 

fl> r (x,t) = f d9e' e/2 (a{6)e ipeX - iEet -a\&)e- ipsX+iE A . (72) 

2 V Jg^Q \ ) 

Then, with G ncss (x,r) = (ip(x,T)ip(0, 0)) Pness and G 1 £ ss (x,t) = (^' r (x, r)^(0, 0)) Pness , the non- 
equilibrium KMS relation for two-point fermion correlation functions is G ness (x, r) = — G l ncss (x, r+ 
Pi) — Gn ess (x,r + fir). This is a non-local relation, hence does not (a priori) have a natural ge- 
ometric interpretation. 

Yet when considering instead the two-point function 

g ness (x,T) := (^(x,T)fi + (0,0)) Pncss 

such a generalized KMS relation is useful as it allows us to establish analytic properties of 
form factors. Indeed, similar arguments, using the exchange relations with the disorder field ji + 
(recall that it has a branch cut towards its right), lead in the non-equilibrium steady state to 
the non-local relations 

g ncss (x, t) = sgn(x) (gi^x, r + fa) + gl css {x, r + /3 r j\ . (73) 

Here we can use a form factor expansion like (|50|) , since the field ip one has only one-particle form 
factors that are nonzero. In particular, in the present case of the non-equilibrium steady state 
(|64|) . integrals exist on the real line, and the large-distance expansion can be obtained exactly. 
The non-equilibrium KMS relation (|73p combined with the form factor expansion then leads 
to nontrivial conditions on the analyticity properties of form factors; in particular, it implies 
the presence of a branch cut, with prescribed jump. As we show in Subsection 15.41 this is in 
agreement with (|65p . 



5 Calculations 

5.1 Mixing 

First, note that the pure-state limit of form factors can be defined by taking the limit V{9) — > oo 
(uniformly on 9) of mixed-state form factors; recall ([25]). By using the cyclic property of the 
trace in order to bring all operators a{9) to the left while keeping all operators a) {9) to the right 
of O, we see that 

lim f£° >+ _(9 1 ,...,9 j ,9 j+1 ,...,9 N ) = (9 N ,...,9 j+1 \0\9 1 ,...,9 J ) (74) 

where there are j plus signs and N — j minus signs as indices on the left-hand side. This holds 
for all values of rapidities such that 9j ^ 9k for j ^ k, and thanks to ([35]) we may evaluate the 
limit for other choices of signs. This suggests that matrix elements of O in % with rapidities 
both on the right and on the left be gathered into the function 

/ew^i.-^M := Jim f e p ;°..,e N (0 1 ,...,9 N ). (75) 
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In particular, the usual form factors, with excited states on the right and the vacuum on the 
left, are f o (0 lt ...,6 N ) = /£...,+ ...,6 N ). 

It will be convenient in order to express formally our result to use a notation for form factors 
with general Liouville states \A) P : 

f p '°[\A) p ] = p (vac\O e \A) p . (76) 

We will also use the notation 

f°[\A) p ] := limJ p '°[\A)P]. (77) 

Proposition 5.1 There is a linear map 

it : End(ft) -> End(ft) 

O U(O) (78) 

satisfying 

U(O t )=H(<9) t > (79) 
and an operator U on C p , such that for all \ A) P G £ p and all O S End('H), 

\U(0)) P = U\0) P , (80) 

and 

r[|A)"]=/^[|^]=/°[ut|^]. (81) 

We have 

U^pf/*^]. (82) 

Further, if O is a local fermionic descendant of the identity (i.e. the identity itself, the fermion 
fields ip(x) and ip(x), their x- derivatives, or normal- ordered products of such), then ii(0) is a 
linear combination of such local fermionic descendants with equal and lower numbers of fermion 
fields. 

In order to prove Proposition 15, 1\ let us consider the operators 

O = a\9 n ) ■ ■ ■ a\9 k+1 )a{e k ) ■ ■ ■ a(6 1 ) =: a^) (83) 

i 

for fixed k, n. These, for all nonnegative integers n > k, form a basis in End(%) according to 
our assumptions. 

Recall the canonical creation and annihilation operators a e (0) and a|(0) on C p , with anti- 
commutation relations (|29p . Let us introduce the normal-ordering operation ° • ° on C p which 
takes the operators a e (0) to the right of the operators a|(0) in the usual way. Then, using (|122|) . 
we have 

:^;|vac>" = n^)l vac ) p ( 84 ) 

i 

"<vac|°<^° = P(vac|n^|^ (85) 
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where 

C e (9) :=l + e* w W. 
Using (j2SJ), (f3"0]) and ([55]) we immediately see that 

"<vac| ° o e ° \e h . . . , OjtOj+u e N ) p + ^ !+! _^_ = (o N , . . . , e j+1 \o\9 u e 3 ) (87) 

(in particular, both sides are nonzero if and only if k = j and n = N). This implies (see (|77p ) 

= p (vac|°(^°|.4> p (88) 

for all A G End (ft). 
From (|82p . we find 

Ua^^U" 1 = 4(9) + e |^a_ e (0). (89) 



Hence, using m 



k+1 1 4(190 , a„(^)\ A f a -(^) a +( 



v c_(eo c+ceoy 11 \c+(ei) c + (e t ) 



|vac) p 



i=n 



c-(9i) c + (9i) j 11 ^c+^i) c_(fli) y 1 ' 

= 0"\vac)P = \oy (90) 

where we used the facts that {at(#), a + (0')} = and that a_|_(#)|vac) p = 0. Taking (f80|) as a 
definition of the map it, we have found 

O l \w&c) p = °ll(0/°|vac) p . (91) 

From (j36[) we observe that ( p (A|il(C>)) p )* = p (i|U(C>)t) p where 1 is obtained by expressing 
\A) P in the Liouville particle basis, and in every term conjugating the coefficient, inverting the 
order of the rapidities, and flipping the particle types (note that this does not give |^) p ). Since 
U is invariant under these operations, we find (P(A\ii(0))P)* = (p(A\XJ\0)p)* = p{A\\J\0^)p = 
P(A\il(0^))P, which implies (|79p. In particular, from (|9ip . we then have 

p (vac\0 £ = P(v&c\iil{0) e i (92) 

Putting together ((55]) and ((92]) we obtain 

f^[\A)P] = p (vac|°il(0) f °|A) P = "(vac|OV) P (93) 



which shows the first equation in (]8ip . The second is obtained similarly but using the hermitian 
conjugate of ((90]) . 

f°[\J^\A) p ] = p (vac|°0 £ °U t |A) p = P(vac\O e \A)P. (94) 

Finally, in order to show the last part of the proposition, we only need to show locality (it 
is clear that U reduces the number of fermion operators, and in fact that it preserves the parity 



28 



of this number). For this, we note that using (|122p and (|89p . given any linear combination 
A = fd8 (u{6)a^{6) + v{6)a{6)), we find 

lT 1 A e V = A t + R, V- l A r V = A r + H{-l) nt - nr (95) 

for some operator R (that depends on u and v). Hence, 

\J- 1 (A e - A r {-l) nl - nT )V = A e - A r {-l) nl - nr '. (96) 

That is, with [j4, •] being a commutator if • is bosonic and an anti-commutator if it is fermionic, 
we deduce 

\[A,il(0)]) p = (A e - A r (-1)"" - nr )\ii(0)) p = U(A e - A r {-l) ni - nr )\0) p = |il([A, O]))" '. (97) 

This gives (|41|). Let us now take = O(x') to be a normal-ordered product of fermion 
fields at the point x' . Then [ip(x),0(x')] = [ifi'(x), 0(x')\ = for all x ^ x'. Hence also 
[tp(x),H(0(x'))} = [^(x), ii(0(x'))} = for all x / x' . Given that il(0(z')) is composed of 
normal-ordered terms each with finitely many factors of creation / annihilation operators, it 
must be a linear combination of normal-ordered product of fermion fields at the point x', plus 
possibly a term proportional to the identity. By translation covariance, this latter term is inde- 
pendent of x 1 . It is clear from the form of U that the term with the highest number of fermion 
field is the original field itself. By dimensional analysis, we have 

il(0(x)) = 0{x) + rn d °- d *F%y{x) (98) 

where are normal-ordered products of fermion fields or the identity field (not involving ex- 
plicitly the mass), are their dimensions and do is the dimension of O; the sum is finite. In 
the limit where V — > oo, the coefficients Fq tend to zero. 



5.2 Mixed-state form factors of twist fields 



We prove Proposition 13.11 by establishing a system of functional differential equations for the 
one- and two-particle form factors of the field ^+ and <r+ (something similar can be done for 
<r_), and by showing that (|42p . (|43p satisfy it. The differential equations are obtained using, 
in particular, the fact that Wick's theorem, as described after Proposition 13. 11 applies for mixed- 
state form factors with higher numbers of particles (and also for traces with insertion of one 
twist field and many creation / annihilation operators). Wick's theorem is a consequence of the 
fact that the twist fields a and ji are simply related to normal-ordered exponentials of bilinear 
expressions in the creation / annihilation operators. For instance the order field <r+ is (a) times 



exp 



J d9 1 d9 2 F ei>e2 {e l ,e 2 )a^{9 x )a^{9 2 ) 



(99) 



where F ei|£2 (#i, #2) = — 2^7 (^l? are essentially the matrix elements on the Hilbert space. 



Hence the equations we derive hold in general for operators of the form (|99p . The differential 
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equations are nonlinear first-order differential equations for the form factors seen as functionals 
of the function W : 9 \— > W{9) that characterizes the density matrix. With the initial condition 
at V{9) = He(W(9)) = oo, given by the matrix elements of the twist fields on the Hilbert space, 
the solution is unique. Indeed, these matrix elements fully characterize the kernel F ei ^ 2 (9i,9 2 ) 
in the bilinear expression. 



5.2.1 Nonlinear functional differential system of equations 

Let (see Q2ZJ) 

m:=Q e (9)^±^M, f ei , e2 (9 1} 9 2 ) := Q^M ^ ^ ^ 

and similarly for higher numbers of insertions of creation and annihilation operators. Using 
5p/6W(P) = -at(fi)a(P)p, we find 

5 Tr (pp + a%9)) _ Tr (p p + a*(9)a^(3)a(f3)) Tr (p p + a%9)) Tr (pa + afi{P)a(/3)) 



8W(/3) Tr(pa + ) Tr(pa+) Tr(pa + ) Tr (pa+) 

Simplifying the right-hand side through Wick's theorem, 

/ e (0)/ + ,_(/3,/3) - fl + ^(0,(3,(3) = / + (/3)/ e ,_(0,/3) - /_(/3)A,+(#,/3), 
this implies, from the definition (|27|h 

S e5(/3-6) \ ~ f + (P)f e ,_(e,P)-f-(l3)f e , + {6,!3) 

+ -, , _ rW (R\ J A") = : — . o w(r\ • l llJ1 ) 



SW{fi) l + e^W) JeK 4cosh 2 ™ 
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Likewise, differentiating Tr (pcr+ a ei (9i)a e2 (9 2 )) /Tr(pcr+), we find 
5 eiStf - gi) e 2 S(?-9 2 ) \ ~ 



5W(f3) l + e £ W) i + e e2W(p) 



4 cosh 2 ™ 



(102) 



Equation (|102|) gives a continuous family of non-linear first-order functional differential equations 
for the W-functionals f ++ (9 1 ,9 2 ), f+-(9 l ,9 2 ), f_+(9 1 ,6 2 ) and f—{9 1 ,9 2 ) (9 U 9 2 G R). Once 
this system of equations is solved, the solution can be fed into (|102p to provide a continuous 
family of linear first-order differential equations for the W-functionals f+{9) and f~(9) for all 
6 € R. The continuous families can be slightly simplified using anti-symmetry under interchange 
of the particles (which is consistent with (|10ip and f|102[) ^ . 

This system of equations can be translated into a system for the mixed-state form factors 
themselves using ([3"7|L since we have, using the notation f e (9) '■= (o")~ 1 /f' M+ (0), f £1 €2 (9\,9 2 ) := 

(a)- 1 f e Zi(e 1 ,9 2 ), 

f e (9) = f € (9), / eijea (0i,0 2 ) = fe U e 2 (9i,0 2 ) + (l + e~^ w ^) 5 ei+€2>0 S(9 1 - 9 2 ). (103) 
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This gives 



sw{p) 



f + (P)f ei -(e,p)-f-(p)f e>+ (o,p) 

4 cosh 2 ™ 



(104) 



6W(J3) 



/ ei , + (g 1 ,^)/ e2 ,_(g 2 ,/3)-/ ei ,_(e 1 ,^)/ e2 , + (g 2 ,/3) 

4 cosh 2 ™ 



(105) 



We see that the delta-function terms have all disappeared. 
5.2.2 Unicity 

From the trace definition of mixed-state form factors, or more explicitly from the expression 
(|81[) of mixed-state form factors in terms of ordinary form factors, we can deduce the form of 
the large-Re(VF(6')) expansion. We have 



where {e} = e or ei,£2, and {8} = 8 or 61,62 • = ei,C2- Then, equations (|104p . (|105|) provide 
uniquely the solutions order by order, once the zeroth order has been fixed. For instance, 



Since the expressions for mixed-state form factors given in Proposition 13. II agree with the usual 
pure-state form factors at Re(W(0)) = 00, then we only have to verify that they satisfy the 
system of equations (|104|) . (|105p in order to prove that they are correct. 

5.2.3 Solution 

In order to verify that the expressions for mixed states form factors of Proposition 13.11 give rise 
to a solution to (|104p and f|105 jl . it is sufficient to take 6\ and 82 in their analyticity region, as 
the system of equations (|lU4j) . (|105p is itself analytic. Hence we consider 9\ S If and 82 G l£, 
see (|4"5]). We then have 



/{«}({*}) = /$({(>})+ J dp ie - w W f { { t\({8},Pi)+ J d{hdP2e- w W- w Wf®({e},Pi,{h)+ 



W) 



XK{6) 





Jw) ht{9) 



2m sinhVF^) sinh(6l - (3) 



e 1 1 



ht \8). 



(106) 



Consider (|104p . On the left-hand side, we find 



e 1 1 1 

2m sinh W(P) sinh(0 - p) 



ht(8) 
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and on-the right-hand side, using (|37|). 

-J- * f ftanh^V 6 - ftanh^Vl 
47risinhW(/3) I V 2 / V 2 / / V* 



These are equal thanks to the relation 



1 / a; a; 
coth tanh 



sinhx 2 V 2 2 
Similarly, consider (|105|) . On the left-hand side, we find 

1 1 / ei e 2 



2vri sinh W{(3) Ysinh(6>i - (3) sinh(6> 2 - P) J \ 
whereas on the right-hand side, again using (1471) 
1 1 



+ -X77T 5x1 f tanh ^T^ 1 * h ei (9i)h e2 (0 2 ) 



Am sinhW(/3) 

x ( tanh P - J ( tanh P - J - ( tanh H - J ( tanh H - J x 

X heMheM. (107) 

Again, these are equal thanks to the relations 

H — , T7 ] tanh - 2 



sinh(#i - /3) sinh(0 2 - /3) / 2 

= — (tanh — coth — — tanh — coth — 

2 V 2 2 2 

i i \ e 2 -e 1 

cotn ■ 



sinh(6»i - p) sinh(6l 2 - 0) J 2 

= — | coth — coth — — tanh — tanh — 

2 V 2 2 2 2 

The latter relations can be verified by looking at both sides as analytic functions of j3. In the 
first relation, both sides have poles at f3 = 9\ and j3 = 9 2 with residues tanh Sl ~® 2 ; and in the 
second the residues are ± coth Sl ~ S2 respectively. In all cases, both sides change sign under 
(3 i — y f3 -\- Z7T and no other poles than those mentioned are found in any strip of width m. 

5.3 A thermodynamic argument 

Consider the one-particle mixed-state form factor 

Tr(p) 
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The idea for the thermodynamic argument is to evaluate the traces by taking the product over 
the one-particle subspaces of all rapidities. This is clear for the denominator since p is diagonal: 

(l + e -^W) _ 1 

Tr(p) "ik^l + e-iW)- 

For the numerator, the argument is less clear, as /x+ is not diagonal. Nevertheless, we may 
imagine at each rapidity that we look at the region far to the right, where there is a cut. Looking 
at this region will be corresponding to a slight imaginary shift in the rapidity, to be implemented 
below. Looking at the other region should give a factor of 1 by the argument presented here. 

Taking the product over the rapidities and taking into account the sign change due to the 



cut, this should give fle'Cl ~~ e 



-W{9'Y 



on the numerator. However, the presence of the particle 



at the rapidity 8 will affect the density of rapidities 8' over which we take the product, as well 
as what we obtain at 8' = 8. Considering a region at a distance L to the right of the origin, we 
expect there to be a factor representing the effect of the branch point, given by 



(vac\a + \8,8')e- iP( >' L oc tanh 



-ip g ,L 



in evaluating the trace on the 8' subspace, when one particle at 8' is present. The dominant 
contribution will be obtained when this factor is equal to 1; then in particular we do get 1 — 
e —W(9') j n evaluating the trace. At 8' = 8 this is not possible as the tanh function is zero - there 
the denominator is simply equal to 1. Hence we must evaluate 



e>+e 



-W(9') 



+ e 



-W(9') 



exp 



^] log ( tanh 

Q<j=Q ^ 



W(8'] 



with a density obtained from 

1 

I " 



2vrzL 



dlog I tanh 



In fact, near to 8 this argument fails, and the fact that we were looking at the far right region 
must be taken into account. This appears to mean that we must replace 8' by 8' — iO (thus 
making e~ ip$/L decaying at large L), and that we must the sum over all values of 6', including 
the value 8. 

Then, the part corresponding to Apg/ gives an infinite contribution. But if the twist field is 
at position x instead of 0, then we must replace L by L — x, and this corresponds to the finite 
change 



exp 

in agreement with (j5l|) . (1551) . 
The other part gives 

i r .,„ b 



2vr B \ 2 



exp 



2W ^v° g ( coth - 



9' + i0 



log tanh 



W(8') 



This is in agreement with (03]). 

Of course, this argument is far from a satisfying derivation, but it does explain the main 
features of (03]). 
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5.4 Analytic structure from non-equilibrium KMS relation 

We re-write the equality (|73p using a form factor expansion like (|50p (and simplifying common 
factors) : 

EC , piepox-eEgT f iep g x-€T(E e +W n oss(8)) 
/ dOe 9 ' 2 — ——h±J0) = sgn(x)V / dOe 9 ' 2 - h+ £ (d). 

In order to establish the equality, we shift the 0-contour by ein/2 for x > and by —ein/2 for 
x < 0, taking care of surrounding the segment from to ±eiir/2. Let us denote by 



Let us consider for instance the case x > 0. Equality will be obtained if the e = + and e = — 
shifted contours cancel each other except at the points where the extra factor e~ trWncss ^ is 
sgn(x) = 1. That is, we require 



e «r/4£. I + y ] + e^'V? I - y ] = (109) 
except for the values of 8 satisfying 

im/3 ri i sinh0 j 

where the functions g r + [9 + ) an d <?l! (0 — -f-) may have poles. One can check that with 
(|43p . these do have poles at the those values and that (|109|) holds elsewhere. Further, equality 
requires that the parts of the contours surrounding the segments emanating from the origin also 
cancel. On the right-hand side, these four terms are 

»0 rml2 



/■u fllT/Z 

/ d0e fl/ ^(0-O)e ipfla - £fl(r+ r) -|- / ^e 9/2 ^(0 + O)e ipea; - £e ( TH 

/0 /■— i7r/2 

d8e e/2 g r L{6 - oyPe^-E e {r+^) + / d8e 9 l 2 gl{6 + o)e ipfla, -- Bfl ( T+A > (110) 
-i7r/2 JO 

and on the left-hand side they are the same except for setting /3 r j to zero. Equality on the 
segment from to z7r/2, and equality on the segment from to —in/ 2, then give the two 
equations (for e = ± respectively) 

_ )e-^ - + 0)e"^ = £(0 - 0) - <?<(0 + 0) 

which imply 



£(6-0) 1-e 



g\{6-0) i- e -ePrEe' 

This gives exactly (j65[) in the case of /il for e = +, and in the case of hf for e = — (and with 
W ness ). 

A similar analysis can be done with x < 0, giving a pole structure in agreement with (|43p 
and showing that fct(f) has no jump through the segment from to in/2, and that h~^(9) has 
no jump through the segment from to —in/2, again in agreement with our results. 
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5.5 General solution as integral-operator kernel 

The nonlinear functional differential equation (|105p presented in Subsection 15.21 is valid for 
general normal-ordered exponential fields of the form (|99p . However, the solution we presented 
is a special solution, for the twist fields cr+. It is not clear a priori that a general solution will 
have the leg- factor structure found. Nevertheless, we may express the general solution of (|105|) 
as an integral-operator kernel constructed out of the functions F eie2 ($x, 62) involved in (|99p . 
We may proceed as follows. Denote : e s : the operator (|99p . with 



S= E / de 1 de 2 F eue2 (e 1 ,02)a ei (0i)a e2 (8 2 ) 

ei,e 2 J 



and consider its two-particle mixed-state form factors, for convenience here taking the two- 
particle state on the left, eij( f 2 (#i, #2|0 eS 0^l vac ) p - We have 

ei / 2 ^i,02|(:e s :)^|vac^ = ei)( f 2 <0i, 2 |U°e s ' ° |vac>» 

= e J 2 {e 1 ,9 2 \Ue § \vacy 

= e J 2 (ei,e 2 \e ljSxj - 1 \YBcr. 

In the first step we used (|9ip . In the second step we use the left-action expressions (|122p and 
the Liouville-space normal ordering in order to simplify the action on the vacuum, defining 

/dQ df) 

(see ([86]) ). In the last step we used U|vac) p = 0. Using ([89]) we then obtain 



p/ "i,e 2 \(: e s :)<|vac)' = e J 2 (0 U # 2 |e D |vac)^ (112) 



where 



On the right-hand side of (|112p we now have a matrix element of a pure exponential, so that we 
may use standard Bogoliubov-transformation techniques. 

One way to implement such techniques is as follows. Consider a general canonical anti- 
commutation algebra with basis bj, b*-, taking discrete indices j = 1,2,... for simplicity, with 
{b*j,bk} = other anti-commutators being zero. Denote by |0) the associated Fock vacuum. 
Denote by V the column vector divided into two blocks, whose elements are 

^=(61,62,...; bl,b* 2 ,...) T , V* = (bl,b*,...; 61,63,...)- 

For general 2-block by 2-block matrix J, we are interested in evaluating the matrix 

<0|We™ \0)_ W:= (Wn W»\ 4) 



(0\e v * JV \0) ' V 
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where on the right-hand side we extracted the block structure. In particular, for the form factor 
problem at hand, we are looking for W± 2 . We evaluate (|114|) by first solving a more general 
problem, writing 



{0\e v * JV \0)W = lim Ti(VV*e 



P— >oo 



V*JV -PN 



N = -V*a z V 
2 



(115) 



where a z acts on the block space. The trace may be evaluated as usual by using its cyclic 
property, moving V along one cycle, and the commutation relations 



{V,V*} = 1, [N,V] = -a z V, [V*JV,V]=MV, M = a x J T a x -J. 



The result is 



1 + e^ Uz e M ) Tr ( VV*e 



D V*JV-PN 



e^e M Ti(e v * JY e-P N 



We take the large-/? limit, keeping only the divergent terms proportional to e 13 on both sides. 
For this, we use e^ <Jz — > e^P\ where Pi is the projector onto the first block. Hence we find 



Pl e M W = P 1& M 



(e M )nW 12 = (e A/ ) 12 => W 12 = ((e M ) u ) (e M ) 



12 



which is the standard expression. 

We now apply this to our form factor problem. Using (|113p and canonically normalized 
operators 

h *M = ^^4(0), b c (o) = -±=^{B) 



we find 



d0id9 2 



Vc+(0i)c+(0 2 



F £1>e2 (e 1 ,e 2 ) (blie^ + etb-eM) (b\ 2 (e 2 ) + e 2 h^ 2 {e 2 , 



(116) 

This gives rise to the matrix of integral operators J by identifying D = V*JV and then to 
M. In the 4 by 4 form taking into account the blocks discussed above as well as the internal 
particle-type e block structure, we obtain 



M 



x —y —y —x 

z x l x l —z 
t _t 



\ 



— z —x —x z 
V x -y -y -x ) 



(117) 



where the integral operators x, y, z have kernels 

F+- (61,62) -F- + {6 2 ,6{\ 



y{6i,6 2 ) 

z(6 u 6 2 ) 



F++{6\,9 2 ) 

f— {Qi,e 2 ) 

y/C+VJC+fa) 
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and #2) = x(02,Oi) is the transposed kernel. We note that M is nilpotent, M 2 = 0, so 

that its exponential can be calculated trivially, 



By direct calculations one then finds 



where e± t 2 = + is on the first row / column and e\p = — is on the second. 

6 Conclusion and discussion 
6.1 Work done 

Generalizing the thermal (or finite-temperature) form factors of |28U33] . we have defined mixed- 
state form factors and evaluated them explicitly in the Ising model, and we have shown how to 
obtain a full large-distance expansion for mixed-state two-point correlation functions of order and 
disorder fields from these form factors. Our method, like that of [28} 133]. differs in an important 
way from other methods more widely used in the literature on massive integrable models: we 
do not "explicitly" perform the trace defining the two-point function with ordinary form factors 
and matrix elements (necessitating a delicate cancelation of divergencies), but rather we define 
mixed-state form factors as simpler traces, evaluated essentially using the cyclic property, and 
show that these traces can be used as building blocks for correlation functions, paralleling the 
ordinary vacuum Kallen-Lehmann expansion. Our basic idea essentially follows from the GNS 
construction of C*-algebras. 

We evaluated mixed-state form factors using a novel technique, showing how to derive for 
them a system of non-linear functional-differential equations. Knowing the vacuum form fac- 
tors, these equations provide uniquely all mixed state form factors up to normalization. These 
techniques appear similar to techniques used in classical integrable models in order to obtain 
bilinear differential equations for tau- functions, which are however usually associated to corre- 
lation functions instead of form factors (see for instance [69]). But we do not know yet if there 
is a full technical equivalence. Our new technique is a departure from the standard techniques 
based on solving a Riemann-Hilbert problem for vacuum or thermal form factors. It is more 
powerful as it does not require any strong analyticity property for the eigenvalues of the density 
matrix. 

We have indicated how to apply our form factor expansions to some physical situations of 
interest: quantum quenches and thermal-flow non-equilibrium steady states. In particular, we 
have found a particular oscillating behaviour in log(mx) in the latter situation, with a frequency 
determined by the temperatures of the asymptotic baths. 




This gives the form factors via the kernel of this operator 



,£(01,02 10 e s :)'|vac>' = VC + (0 1 )C + (9 2 )Z eue2 (9 1 ,9 2 ) 
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Our form factor expansion can be trivially re-expressed as the determinant of a matrix linear 
integral operator, generalizing for instance the ideas of [70] (as was done in [28] for the thermal 
case). Prom this perspective, we note that techniques of [71] for obtaining such representations 
of thermal correlation functions in the XY model are likely to be applicable as well to the case 
of general diagonal density matrices. 

6.2 Further developments 

There are many directions that the current work has left to be explored. 

First, a development of both applications to quantum quenches and to the non-equilibrium 
energy-flow steady state are clearly desirable. In particular, this will require a full extension of 
our formalism to more general functions W(9), including with discontinuities and with regions 
of negativity. 

Second, and related to the above, it would be most interesting to understand what generalizes 
the "quantization on the circle" viewpoint of thermal correlation functions, and the Matsubara 
frequencies. From our analysis of the quantum quench and non-equilibrium steady state appli- 
cations, it appears that the generalization will be based on a study of the singularity structure of 
the filling factors associated to the density matrix (bosonic or fermionic, depending on the semi- 
locality structure of the field considered). That is, there should be a correspondence between 
the spectrum in an alternative quantization scheme, where mixed-state correlation functions are 
reproduced by vacuum expectation values, and the singularity structure of these filling factors. 
It would also be interesting to see how these principles, here understood in the context of the 
Ising model, can be applied to more general QFT in mixed states. 

Third, by the symmetry arguments of [26] (or generalizing the arguments of |70j), one can 
easily see that the mixed-state correlation functions we considered satisfy a system of bilinear 
partial differential equations which are essentially the Hirota form of the sinh-Gordon equation. 
As was discovered in [33], the thermal form factors of the Ising model can be seen as the initial 
scattering data for the associated inverse scattering problem. We expect the same to hold in 
the cases of more general mixed states. In particular, it is possible that the choices of density 
matrices be in correspondence with the possible initial conditions for the sinh-Gordon equation. 

Fourth, it would be most interesting to analyze the large-time expansion of correlation func- 
tions. The aforementioned sinh-Gordon equation method may be of use, as well as, for instance, 
the virial expansion method of |30| . 

Fifth, we note that although we have not evaluated the normalization of form factors (and in 
the non-equilibrium steady-state case, some subtleties arise), it is likely that our techniques of 
deriving non-linear functional differential equations can be used similarly for the normalization, 
along with a careful cancellation of divergencies, in order to evaluate it exactly. 

Finally, not only our techniques can be immediately applied to any other free-fermion (or free- 
boson) model, like the sine-Gordon model at the free-fermion point, but also we expect similar 
techniques to provide an interesting new approach to the study of thermal and more generally 
mixed-state correlation functions in integrable models with non-trivial scattering matrix. For 
instance, it is likely that one can obtain in a very efficient way the low-temperature expansion 
of thermal correlation functions, and more generally there is a hope, using thermodynamic 
arguments combined with exact functional-differential equations, to obtain exact mixed-state 
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form factors. 
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A Simple properties in the Liouville space 
A.l Left- and right-actions 

The Liouville left-action is defined in (|32p . Similarly, to every B £ End% one can define a 
right-action B r by B r \A) p = \AB) P . The right-action linear map B i— >• B r is an algebra anti- 
homomorphism (AB) r = B r A r . Using Definition (|22p . one can see that on conjugate vectors we 
have 

p (A\B e = P (B^A\, p (A\B r = p (ApB^p- 1 \. (118) 
Clearly, left- and right-action Liouville operators commute with each other, A e, B r = B r A e . 

A. 2 Hermitian conjugations 

It is a simple matter to translate more generally the Hermitian conjugation of operators on H 
onto that of operators on C p (as we mentioned, we denote the Hermitian conjugation in both 
cases by '). In particular, conjugating the equations (|118p . we find, for every B £ End('H), 

(B e y = (B^) e , (ir)t = (pi?y- i ) r . ( ii9) 

That is, the Hermitian conjugation commutes with the left-action map, but not with the right- 
action map. Specializing to operators a e (9) and using 

pa e (9)p- 1 = e- eW{e) a e {6) (120) 

as well as linearity of the right-action map, we obtain 

(a e (0) r ) t = e eW ^a-%e) r . (121) 

A. 3 Creation and annihilation operators 

We can of course express the left- and right-actions of the Hilbert space creation and annihilation 
operators a e {6) in terms of the operators a e (6»), aj(0) ([29]) • Let (-l) n = (-!)/<*? at (<?) a (0) e 
End (7^) be the operator that gives the parity of the number of Hilbert-space particles. Then 
(— l) n<? ~ nr £ End(£p) is the operator that gives the parity of the number of Liouville particles, 
which can also be written (-l) n = (-l)E e /^a|(0)a e (e)/(i+e-^( e ))_ From ^ and it i s c l ea r 

that we should have a^O) 1 = 1+ ^% e) + &a_ £ (0) and a*(8) r (-l f-^ = 1+ *1 ( 1 (9) + ? e a_ e (0), 
where in both cases, the second term must be present in order to account for the fact that no 
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delta-function contact term is obtained in permuting Liouville particles, but that such contact 
terms appear in commuting Hilbert-space annihilation and creation operators. The values of 
q t and q e can be fixed uniquely using the Hermitian structure of C p , more precisely (|119p and 
(|12ip (which imply in particular (n £ )^ = n , (n r )t = n r ). We find 

„ em t a|(g) a_ £ (0) 



l + e -eW(8) 1 + e eW(9) 

a e (9Y(-lf- nr = \ w{e) (at(g) + a_ e (g)). (122) 



1 + e -^(«) 

A. 4 Hamiltonian and momentum operators 



Finally, we note that our choice (|23|) of density matrix p guarantees that space and time trans- 
lation invariance hold as well in C p (however, of course, we do not have in general Poincare 
invariance). Indeed, we can define operators 



J 1 + e -tW(6) 

P != p,_ P , = E /^ £w J^KW (123) 



which commute with each other, are Hermitian, and are diagonalized on the basis . . . , 0zv)ei,...,ejsn 
with in particular H|vac) p = P|vac) p = 0. For any A G End(£p), space and time translations 
are given by 

A(x, t) = e 1 ™- 1 ' 9 * Ae- im+iPx . (124) 

This is in agreement with the left- and right-action maps: with A(x, t) = e tHt ~ lPx Ae~ lHt+lPx 
we have 

A{x,tY = A\x,t), A{x,t) r = A r {x,t) (125) 

where we used the homomorphism and anti-homomorphism properties of the left- and right- 
action maps, respectively, as well as the fact that left- and right-action Liouville operators 
commute with each other. It is also in agreement with the correspondence between Hilbert 
space operators and Liouville space vectors: 

\A{x,t)) p = e im - iPx \A) p . (126) 



B Reality of correlation function 

The order and disorder two-point correlation functions should be real functions. This is a non- 
trivial fact that encompasses information both about the hermitian conjugation properties of 
the fields, and about their locality. Reality is not explicit in the form factor expansion (|59p . We 
provide a simple check here that the one-particle term is indeed real. The verification to higher 
particles is more involved and we defer it to another paper. 
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In order to verify that the one-particle order of the form factor expansion is real, we note 
that 

K t {6f = e eW{e) K_ € {6). 

Hence, after changing variable to ej t— > —ej, the imaginary part of the one-particle integral on 
right-hand side of (|59p is proportional to 

I d6 (l - e- eW{9 ^ K e {6)e iepB ^- ieEst . (127) 

The function J £ (0) = (1- eT eW W)K e (0) = i € ~ l g{6 + ieO)~ 2€ /(2tt) can be analytically continued, 
and satisfies 

J + (9 + ivr/2) = -J_(0 - in/2). 
Hence, shifting the ^-contour by ien/2 (no pole is crossed), we find that (|127p is zero. 

C Non-equilibrium steady-state form factor expansion 



We start with the expression (|57p . In order to obtain convergent integrals, we shift towards 
imaginary directions. Let us consider one term in the sum £^ 6 . It is convenient to shift 
all 0,,-contours with €j = + by a quantity +iir, and not to shift at all the integral contours with 
€k = — . Using crossing symmetry (|49p (more precisely, similar identities for all rapiditiez 9j), 
we see that on the shifted contours the integrand in (|57p becomes 



'N 



'Ef=i(iPe j x-iE gj tj 



where the number of minus signs J is the number of shifted contours (the number of indices j 
with ej = +). Further, the shifting gives rise to residue contributions coming from the poles of 
the factors (1 — e -M/ncss ( 61 -')) -1 , and to an integral running on both sides of the imaginary segment 
between and m because the point 6j = is not analytic (discontinuity of W ncss (9j)). This 
means that we may replace every ^-integral by 



dOj = -) 

— J dOj + residues + imaginary segment (ej = +). 



Hence, under the sum over ex,... , ejv, only the residues and imaginary segment contributions 
remain. By exchanging particles, we may assume that we take residues for the first P particles 
and we integrate on imaginary segments for the remaining Q = N — P particles, taking care of 
putting the extra combinatorial factor N\/(P\Q\). 

The pole in 9j whose residue is taken is at position iir/2 + a(n) for all n G Z \ {0} (see 
(|69|) ). Evaluating the residues and putting together the imaginary segment contributions, after 
a change to an integration over a real variable, one obtains ()67p . 
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In order to evaluate the leading large-ma; behavior (setting t = for simplicity) for the 
disorder two-point function, we consider N = and M = 1, and only the part of the integral 
near and tt is sufficient. Hence we consider 

sinh ( — cos 



/_ \ rTT sinn i o cos c i 

\»/Pno S3 e -xg ncss / dde -mx S in6 ^ ^/i»+(i0)/lL(*0). (128) 



sinh ^ cos #J sinh cc 

The large-mx behavior is evaluated by expanding for 9 near to and tt. We can then use 
h$ + (iO) h^_(i6) oc (i#) 2 * 7 near to 9 = from (|66p . and similarly oc (z7r — i9)~ 2l ~ < near to = tt, 
using (jl8|) . Omitting the overall finite, real (temperature-dependent) factor, we then obtain, 
asymptotically, 

f°° 1 

e iB / d9 e - mx6 e 2 ^ + c.c. oc cos(2 7 log(mx) + B) 

Jo mx 

JB 



for some phase e . A more careful calculation gives (|71|) with 



sinh 

A = 2 s 2 . lr(l + 2t7)| (129) 

sinn Rj— smn i - 2 2 — 

5 = arg(r(l + 2i 7 )) + 

+ 27 j( e|>1 ^ iinM l0g C ° th ~ 2~ + 27 | K1 ^ (sinhtf " g * 

Higher values of M are likely to lead to contributions of exactly the same form. Indeed, we 
will have integrals of the type 

£ d$idB 2 ^tai / 1 e- mx ( sin9l+sine2 ) x leg factors. 

For 9\ ~ and #2 ~ and vice versa, the integrand has a second-order pole. This leads to 
logarithmic divergences, which we expect can be re-absorbed into the normalization of the field. 
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